From 14fb78465fe33598affeb1abd1bfdacd59f16f26 Mon Sep 17 00:00:00 2001 From: masera Date: Mon, 15 Sep 2014 16:38:41 +0200 Subject: [PATCH] Latest version of the CA tracker. Big improvement in terms of CPU speed. --- ITS/UPGRADE/AliITSUAux.h | 6 +- ITS/UPGRADE/AliITSUCACell.h | 106 ++ ITS/UPGRADE/AliITSUCATracker.cxx | 1274 ++++++++++++++++++++++ ITS/UPGRADE/AliITSUCATracker.h | 144 +++ ITS/UPGRADE/AliITSUCATrackingStation.cxx | 368 +++++++ ITS/UPGRADE/AliITSUCATrackingStation.h | 159 +++ ITS/UPGRADE/AliITSUTrackerSA.cxx | 885 --------------- ITS/UPGRADE/AliITSUTrackerSA.h | 111 -- ITS/UPGRADE/AliITSUTrackerSAaux.h | 268 ----- ITS/UPGRADE/CMakelibITSUpgradeRec.pkg | 5 +- ITS/UPGRADE/ITSUpgradeRecLinkDef.h | 2 +- ITS/UPGRADE/macros/testTrackerCA.C | 97 +- 12 files changed, 2149 insertions(+), 1276 deletions(-) create mode 100644 ITS/UPGRADE/AliITSUCACell.h create mode 100644 ITS/UPGRADE/AliITSUCATracker.cxx create mode 100644 ITS/UPGRADE/AliITSUCATracker.h create mode 100644 ITS/UPGRADE/AliITSUCATrackingStation.cxx create mode 100644 ITS/UPGRADE/AliITSUCATrackingStation.h delete mode 100644 ITS/UPGRADE/AliITSUTrackerSA.cxx delete mode 100644 ITS/UPGRADE/AliITSUTrackerSA.h delete mode 100644 ITS/UPGRADE/AliITSUTrackerSAaux.h diff --git a/ITS/UPGRADE/AliITSUAux.h b/ITS/UPGRADE/AliITSUAux.h index abad138c368..59331da3df6 100644 --- a/ITS/UPGRADE/AliITSUAux.h +++ b/ITS/UPGRADE/AliITSUAux.h @@ -21,7 +21,8 @@ using namespace TMath; namespace AliITSUAux { - void BringTo02Pi(double &phi); + template + void BringTo02Pi(F &phi); Bool_t OKforPhiMin(double phiMin,double phi); Bool_t OKforPhiMax(double phiMax,double phi); Double_t MeanPhiSmall(double phi0, double phi1); @@ -44,7 +45,8 @@ namespace AliITSUAux { } //_________________________________________________________________________________ -inline void AliITSUAux::BringTo02Pi(double &phi) { +template +inline void AliITSUAux::BringTo02Pi(F &phi) { // bring phi to 0-2pi range if (phi<0) phi+=TwoPi(); else if (phi>TwoPi()) phi-=TwoPi(); } diff --git a/ITS/UPGRADE/AliITSUCACell.h b/ITS/UPGRADE/AliITSUCACell.h new file mode 100644 index 00000000000..64f30a1ff05 --- /dev/null +++ b/ITS/UPGRADE/AliITSUCACell.h @@ -0,0 +1,106 @@ +#ifndef ALIITSUCACELL_H +#define ALIITSUCACELL_H + +#include +#include + +class AliITSUCACell { +public: + AliITSUCACell(int xx = 0u,int yy = 0u, int zz = 0u, int dd0 = 0, + int dd1 = 0, float curv = 0.f, float n[3] = 0x0) + : f1OverR(curv), fd0(dd0), fd1(dd1), fN(), fVector() + { + fVector.reserve(4); + fVector.push_back(xx); + fVector.push_back(yy); + fVector.push_back(zz); + fVector.push_back(1u); + if(n) { + fN[0] = n[0]; + fN[1] = n[1]; + fN[2] = n[2]; + } + } + + int x() const { return fVector[0]; } + int y() const { return fVector[1]; } + int z() const { return fVector[2]; } + int d0() const { return fd0; } + int d1() const { return fd1; } + int GetLevel() const { return fVector[3]; } + float GetCurvature() const { return f1OverR; } + float* GetN() { return fN; } + + void SetLevel(int lev) { fVector[3] = lev; } + + int operator()(const int i) { return fVector[4 + i]; } + + size_t NumberOfNeighbours() { return (fVector.size() - 4u); } + + bool Combine(AliITSUCACell &neigh, int idd) + { + // From outside inward + if (this->y() == neigh.z() && this->x() == neigh.y()) // Cells sharing two points + { + fVector.push_back(idd); + if (neigh.GetLevel() + 1 > GetLevel()) + { + SetLevel(neigh.GetLevel() + 1u); + } + return true; + } + return false; + } + +private: + float f1OverR; + int fd0,fd1; + float fN[3]; + std::vector fVector; +}; + +class AliITSUCARoad { +public: + AliITSUCARoad() : Elements(), N(0) + { + ResetElements(); + } + + AliITSUCARoad(int layer, int idd) : Elements(), N() + { + ResetElements(); + N = 1; + Elements[layer] = idd; + } + + AliITSUCARoad(const AliITSUCARoad& copy) : Elements(), N(copy.N) + + { + for ( int i=0; i<5; ++i ) + { + Elements[i] = copy.Elements[i]; + } + } + + int &operator[] (const int &i) { + return Elements[i]; + } + + void ResetElements() + { + for ( int i=0; i<5; ++i ) + Elements[i] = -1; + N = 0; + } + + void AddElement(int i, int el) + { + ++N; + Elements[i] = el; + } + + int Elements[5]; + int N; +}; + +#endif // ALIITSUCACELL_H diff --git a/ITS/UPGRADE/AliITSUCATracker.cxx b/ITS/UPGRADE/AliITSUCATracker.cxx new file mode 100644 index 00000000000..18faadb5e34 --- /dev/null +++ b/ITS/UPGRADE/AliITSUCATracker.cxx @@ -0,0 +1,1274 @@ +//-*- Mode: C++ -*- +// ************************************************************************** +// This file is property of and copyright by the ALICE ITSU Project * +// ALICE Experiment at CERN, All rights reserved. * +// * +// Primary Author: Maximiliano Puccio * +// for the ITS Upgrade project * +// * +// Permission to use, copy, modify and distribute this software and its * +// documentation strictly for non-commercial purposes is hereby granted * +// without fee, provided that the above copyright notice appears in all * +// copies and that both the copyright notice and this permission notice * +// appear in the supporting documentation. The authors make no claims * +// about the suitability of this software for any purpose. It is * +// provided "as is" without express or implied warranty. * +// * +//*************************************************************************** + +#include +#include +#include +#include +#include +#include "AliLog.h" +#include "AliESDEvent.h" +#include "AliITSUClusterPix.h" +#include "AliITSUCATracker.h" +#include "AliITSUReconstructor.h" +#include "AliITSURecoDet.h" +#include "AliESDtrack.h" +#include + +using TMath::Abs; +using TMath::Sort; +using TMath::Sqrt; +using std::sort; +using std::cout; +using std::endl; +using std::flush; + +ClassImp(AliITSUCATracker) + + +// tolerance for layer on-surface check +const Double_t AliITSUCATracker::fgkToler = 1e-6; +const Double_t AliITSUCATracker::fgkChi2Cut = 600.f; +const int AliITSUCATracker::fgkNumberOfIterations = 1; +const float AliITSUCATracker::fgkR[7] = {2.33959,3.14076,3.91924,19.6213,24.5597,34.388,39.3329}; +// +const float kmaxDCAxy[5] = {0.05f,0.04f,0.05f,0.2f,0.4f}; +const float kmaxDCAz[5] = {0.2f,0.4f,0.5f,0.6f,3.f}; +const float kmaxDN[4] = {0.002f,0.009f,0.002f,0.005f}; +const float kmaxDP[4] = {0.008f,0.0025f,0.003f,0.0035f}; +const float kmaxDZ[6] = {0.1f,0.1f,0.3f,0.3f,0.3f,0.3f}; +//const float kSigma2 = 0.0005f * 0.0005f; + +//__________________________________________________________________________________________________ +static inline float invsqrt(float _x ) +{ + // + // The function calculates fast inverse sqrt. Google for 0x5f3759df. + // Credits to ALICE HLT Project + // + + union { float f; int i; } x = { _x }; + const float xhalf = 0.5f * x.f; + x.i = 0x5f3759df - ( x.i >> 1 ); + x.f = x.f * ( 1.5f - xhalf * x.f * x.f ); + return x.f; +} + +//__________________________________________________________________________________________________ +static inline float Curvature(float x1, float y1, float x2, float y2, float x3, float y3) +{ + // + // Initial approximation of the track curvature + // +// return - 2.f * ((x2 - x1) * (y3 - y2) - (x3 - x2) * (y2 - y1)) +// * invsqrt(((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1)) * +// ((x2 - x3) * (x2 - x3) + (y2 - y3) * (y2 - y3)) * +// ((x3 - x1) * (x3 - x1) + (y3 - y1) * (y3 - y1))); + + //calculates the curvature of track + float den = (x3 - x1) * (y2 - y1) - (x2 - x1) * (y3 - y1); + if(den * den < 1e-32) return 0.f; + float a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den; + if((y2 - y1) * (y2 - y1) < 1e-32) return 0.f; + float b = -(x2 * x2 - x1 * x1 + y2 * y2 - y1 * y1 + a * (x2 - x1)) / (y2 - y1); + float c = -x1 * x1 - y1 * y1 - a * x1 - b * y1; + float xc = -a / 2.f; + + if((a * a + b * b - 4 * c) < 0) return 0.f; + float rad = TMath::Sqrt(a * a + b * b - 4 * c) / 2.f; + if(rad * rad < 1e-32) return 1e16; + + if((x1 > 0.f && y1 > 0.f && x1 < xc)) rad *= -1.f; + if((x1 < 0.f && y1 > 0.f && x1 < xc)) rad *= -1.f; + // if((x1<0 && y1<0 && x10 && y1<0 && x1 kAlmost0 ? +// 0.5f * (k1 * k2 * (y1 - y3) + k2 * (x1 + x2) - k1 * (x2 + x3)) / (k2 - k1) : 1e12f; +//} + +//__________________________________________________________________________________________________ +static inline bool CompareAngles(float alpha, float beta, float tolerance) +{ + const float delta = TMath::Abs(alpha - beta); + return (delta < tolerance || TMath::Abs(delta - TMath::TwoPi()) < tolerance); +} + +//__________________________________________________________________________________________________ +AliITSUCATracker::AliITSUCATracker(AliITSUReconstructor* rec) : + fReconstructor(rec), + fITS(0), + fMatLUT(0), + fUseMatLUT(kFALSE), + fCurrMass(0.14), + fLayer(), + fUsedClusters(), + fChi2Cut(fgkChi2Cut), + fPhiCut(1), + fZCut(0.5f), + fCandidates() +{ + // This default constructor needs to be provided + if (rec) Init(rec); + for (int i = 0; i < 4; ++i) + { + fCandidates[i] = new TClonesArray("AliITSUTrackCooked",100000); + } + +#ifdef _TUNING_ + for (int i = 0; i < 6; ++i) + { + fGDZ[i] = new TH1F(Form("DGZ%i",i),Form("DZ%i;#Deltaz;Entries",i),500,0.f,5.f); + fGDXY[i] = new TH1F(Form("DGPhi%i",i),Form("#Delta#Phi%i;#DeltaPhi;Entries",i),500,0.f,TMath::TwoPi()); + fFDZ[i] = new TH1F(Form("DFZ%i",i),Form("DZ%i;#Deltaz;Entries",i),500,0.f,5.f); + fFDXY[i] = new TH1F(Form("DFPhi%i",i),Form("#Delta#Phi%i;#Delta#Phi;Entries",i),500,0.f,TMath::TwoPi()); + } + for (int i = 0; i < 5; ++i) + { + fGDCAZ[i] = new TH1F(Form("DCAGZ%i",i),Form("DCAZ%i;#Deltaz;Entries",i),500,0.f,5.f); + fGDCAXY[i] = new TH1F(Form("DCAGXY%i",i),Form("DCAXY%i;#Deltar;Entries",i),500,0.f,5.f); + fFDCAZ[i] = new TH1F(Form("DCAFZ%i",i),Form("DCAZ%i;#Deltaz;Entries",i),500,0.f,5.f); + fFDCAXY[i] = new TH1F(Form("DCAFXY%i",i),Form("DCAXY%i;#Deltar;Entries",i),500,0.f,5.f); + } + for(int i = 0; i < 4; ++i) + { + fGoodCombChi2[i] = new TH1F(Form("goodcombchi2%i",i),Form("%i;#chi^{2};Entries",i),2000,0,0.1); + fFakeCombChi2[i] = new TH1F(Form("fakecombchi2%i",i),Form("%i;#chi^{2};Entries",i),2000,0,0.1); + fGoodCombN[i] = new TH1F(Form("goodcombn%i",i),Form("%i;#Deltan;Entries",i),300,0.f,0.03f); + fFakeCombN[i] = new TH1F(Form("fakecombn%i",i),Form("%i;#Deltan;Entries",i),300,0.f,0.03f); + } + fGoodCombChi2[4] = new TH1F("goodcomb4",";#chi^{2};Entries",200,0,500); + fFakeCombChi2[4] = new TH1F("fakecomb4",";#chi^{2};Entries",200,0,500); + fTan = new TH1F("tan","tan",2500,0,0.5); + fTanF = new TH1F("tanF","tanF",2500,0,0.5); + fPhi = new TH1F("phi","phi",2500,0,TMath::Pi()); + fPhiF = new TH1F("phi","phiF",2500,0,TMath::Pi()); + fNEntries = new TH1F("nentries","nentries",2001,-0.5,2000.5); +#endif +} + +//__________________________________________________________________________________________________ +AliITSUCATracker::AliITSUCATracker(const AliITSUCATracker &t): AliTracker(t), + fReconstructor(t.fReconstructor), + fITS(t.fITS), + fMatLUT(t.fMatLUT), + fUseMatLUT(t.fUseMatLUT), + fCurrMass(t.fCurrMass), + fLayer(), + fUsedClusters(), + fChi2Cut(fgkChi2Cut), + fPhiCut(), + fZCut(0.5f), + fCandidates() +{ + // The copy constructor is protected +} + +//__________________________________________________________________________________________________ +AliITSUCATracker::~AliITSUCATracker() +{ + // d-tor + for (int i = 0; i < 4; ++i) + { + if (fCandidates[i]) + delete fCandidates[i]; + } + delete fMatLUT; +} + +//__________________________________________________________________________________________________ +void AliITSUCATracker::Init(AliITSUReconstructor* rec) +{ + // init with external reconstructor + // + fITS = rec->GetITSInterface(); + // + // create material lookup table + const int kNTest = 1000; + const double kStepsPerCM = 5; + fMatLUT = new AliITSUMatLUT(fITS->GetRMin(),fITS->GetRMax(), \ + Nint(kStepsPerCM * (fITS->GetRMax() - fITS->GetRMin()))); + double zmn = 1e6; + for (int ilr = fITS->GetNLayers(); ilr--;) + { + AliITSURecoLayer* lr = fITS->GetLayer(ilr); + if (zmn>Abs(lr->GetZMin())) zmn = Abs(lr->GetZMin()); + if (zmn>Abs(lr->GetZMax())) zmn = Abs(lr->GetZMax()); + } + fMatLUT->FillData(kNTest,-zmn,zmn); + // +} + +//__________________________________________________________________________________________________ +void AliITSUCATracker::CellsTreeTraversal(vector &roads, + const int &iD, const int &doubl) +{ + + if (doubl < 0) return; + + roads.back().AddElement(doubl,iD); + //Road tmp = roads.back(); + const int currentN = roads.back().N; + for (size_t iN = 0; iN < fCells[doubl][iD].NumberOfNeighbours(); ++iN) + { + const int currD = doubl - 1; + const int neigh = fCells[doubl][iD](iN); + + if (iN > 0) + { + roads.push_back(roads.back()); + roads.back().N = currentN; + } + + CellsTreeTraversal(roads,neigh,currD); + } + + fCells[doubl][iD].SetLevel(0u); // Level = -1 +} + +//__________________________________________________________________________________________________ +Int_t AliITSUCATracker::Clusters2Tracks(AliESDEvent *event) +{ + // This is the main tracking function + // The clusters must already be loaded + + int ntrk = 0, ngood = 0; + for (int iteration = 0; iteration < fgkNumberOfIterations; ++iteration) + { + + fCandidates[0]->Clear(); + fCandidates[1]->Clear(); + fCandidates[2]->Clear(); + fCandidates[3]->Clear(); + + FindTracksCA(iteration); + + for (int iL = 3; iL >= 0; --iL) + { + const int nCand = fCandidates[iL]->GetEntries(); + int index[nCand]; + float chi2[nCand]; + + for (int iC = 0; iC < nCand; ++iC) + { + AliITSUTrackCooked *tr = (AliITSUTrackCooked*)fCandidates[iL]->At(iC); +// Int_t clInfo[2 * AliITSUAux::kMaxLayers]; +// for (unsigned int i = 0; i < 2 * AliITSUAux::kMaxLayers; ++i) +// { +// clInfo[i] = -1; +// } +// for (int k = 0; k < tr->GetNumberOfClusters(); ++k) +// { +// const int layer = (tr->GetClusterIndex(k) & 0xf0000000) >> 28; +// const int idx = (tr->GetClusterIndex(k) & 0x0fffffff); +// clInfo[layer << 1] = idx; +// } + chi2[iC] = tr->GetChi2();//(RefitTrack(tr,clInfo,0.,-1)); + index[iC] = iC; + CookLabel(tr,0.f); + } + + TMath::Sort(nCand,chi2,index,false); + + for (int iUC = 0; iUC < nCand; ++iUC) + { + const int iC = index[iUC]; + if (chi2[iC] < 0.f) + { + continue; + } +// if (chi2[iC] > fChi2Cut) +// { +// break; +// } + + AliITSUTrackCooked *tr = (AliITSUTrackCooked*)fCandidates[iL]->At(iC); + + bool sharingCluster = false; + for (int k = 0; k < tr->GetNumberOfClusters(); ++k) + { + const int layer = (tr->GetClusterIndex(k) & 0xf0000000) >> 28; + const int idx = (tr->GetClusterIndex(k) & 0x0fffffff); + if (fUsedClusters[layer][idx]) + { + sharingCluster = true; + break; + } + } + + if (sharingCluster) + continue; + + for (int k = 0; k < tr->GetNumberOfClusters(); ++k) + { + const int layer = (tr->GetClusterIndex(k) & 0xf0000000) >> 28; + const int idx = (tr->GetClusterIndex(k) & 0x0fffffff); + fUsedClusters[layer][idx] = true; + } + + AliESDtrack outTrack; + CookLabel(tr,0.f); + ntrk++; + if(tr->GetLabel() >= 0) + { + ngood++; +#ifdef _TUNING_ + fGoodCombChi2[4]->Fill(chi2[iC]); + } + else + { + fFakeCombChi2[4]->Fill(chi2[iC]); +#endif + } + + outTrack.UpdateTrackParams(tr,AliESDtrack::kITSin); + outTrack.SetLabel(tr->GetLabel()); + event->AddTrack(&outTrack); + } + } + } + Info("Clusters2Tracks","Reconstructed tracks: %d",ntrk); + if (ntrk) + Info("Clusters2Tracks","Good tracks/reconstructed: %f",Float_t(ngood)/ntrk); +// + return 0; +} + +//__________________________________________________________________________________________________ +void AliITSUCATracker::FindTracksCA(int) +{ + +#ifdef _TUNING_ + unsigned int numberOfGoodDoublets = 0, totalNumberOfDoublets = 0; + unsigned int numberOfGoodCells = 0, totalNumberOfCells = 0; + unsigned int cellsCombiningSuccesses = 0, cellsWrongCombinations = 0; + unsigned int totalNumberOfRoads = 0; +#endif + + vector dLUT[5]; + for (int iL = 0; iL < 6; ++iL) { + if (iL < 5) { + dLUT[iL].resize(fLayer[iL + 1].GetNClusters(),0); + } + for (int iC = 0; iC < fLayer[iL].GetNClusters(); ++iC) { + ClsInfo_t* cls = fLayer[iL].GetClusterInfo(iC); + const float tanL = (cls->z - GetZ()) / cls->r; + const float extz = tanL * (fgkR[iL + 1] - cls->r) + cls->z; + const int nClust = fLayer[iL + 1].SelectClusters(extz - 5 * fZCut, extz + 5 * fZCut, + cls->phi - fPhiCut, cls->phi + fPhiCut); + bool first = true; + + for (int iC2 = 0; iC2 < nClust; ++iC2) { + const int iD2 = fLayer[iL + 1].GetNextClusterInfoID(); + ClsInfo_t* cls2 = fLayer[iL + 1].GetClusterInfo(iD2); + const float dz = tanL * (cls2->r - cls->r) + cls->z - cls2->z; + if (TMath::Abs(dz) < kmaxDZ[iL] && CompareAngles(cls->phi, cls2->phi, fPhiCut)) { + if (first && iL > 0) { + dLUT[iL - 1][iC] = fDoublets[iL].size(); + first = false; + } + const float dTanL = (cls->z - cls2->z) / (cls->r - cls2->r); + const float phi = TMath::ATan2(cls->y - cls2->y, cls->x - cls2->x); + fDoublets[iL].push_back(Doublets(iC,iD2,dTanL,phi)); +#ifdef _TUNING_ + if (fLayer[iL].GetClusterSorted(iC)->GetLabel(0) == fLayer[iL + 1].GetClusterSorted(iD2)->GetLabel(0) && + fLayer[iL].GetClusterSorted(iC)->GetLabel(0) > 0) { + numberOfGoodDoublets++; + fGDZ[iL]->Fill(dz); + fGDXY[iL]->Fill(fabs(cls->phi - cls2->phi)); + } else { + fFDZ[iL]->Fill(dz); + fFDXY[iL]->Fill(fabs(cls->phi - cls2->phi)); + } + totalNumberOfDoublets++; +#endif + } + } + fLayer[iL + 1].ResetFoundIterator(); + } + } + + vector tLUT[4]; + for (int iD = 0; iD < 5; ++iD) + { + if (iD < 4) { + tLUT[iD].resize(fDoublets[iD + 1].size(),0); + } + for (size_t iD0 = 0; iD0 < fDoublets[iD].size(); ++iD0) + { + const int idx = fDoublets[iD][iD0].y; + bool first = true; + for (size_t iD1 = dLUT[iD][idx]; idx == fDoublets[iD + 1][iD1].x;++iD1) + { + //cout << dLUT[iD][fDoublets[iD][iD0].y] << " " << dLUT[iD][fDoublets[iD][iD0].y + 1] << " " << fLayer[iD + 1].GetNClusters() << endl; + if (TMath::Abs(fDoublets[iD][iD0].tanL - fDoublets[iD + 1][iD1].tanL) < 0.025 && // TODO: cuts as parameters + TMath::Abs(fDoublets[iD][iD0].phi - fDoublets[iD + 1][iD1].phi) < 0.14) { // TODO: cuts as parameters + const float tan = 0.5f * (fDoublets[iD][iD0].tanL + fDoublets[iD + 1][iD1].tanL); + const float extz = -tan * fLayer[iD][fDoublets[iD][iD0].x]->r + fLayer[iD][fDoublets[iD][iD0].x]->z; + if (fabs(extz - GetZ()) < kmaxDCAz[iD]) { +#ifdef _TUNING_ + fGood = (fLayer[iD].GetClusterSorted(fDoublets[iD][iD0].x)->GetLabel(0) == fLayer[iD + 1].GetClusterSorted(fDoublets[iD][iD0].y)->GetLabel(0) && + fLayer[iD].GetClusterSorted(fDoublets[iD][iD0].x)->GetLabel(0) == fLayer[iD + 2].GetClusterSorted(fDoublets[iD + 1][iD1].y)->GetLabel(0) && + fLayer[iD].GetClusterSorted(fDoublets[iD][iD0].x)->GetLabel(0) > 0); +#endif + float curv, n[3]; + if (CellParams(iD, fLayer[iD][fDoublets[iD][iD0].x], fLayer[iD + 1][fDoublets[iD][iD0].y], + fLayer[iD + 2][fDoublets[iD + 1][iD1].y], curv, n)) { + if (first && iD > 0) { + tLUT[iD - 1][iD0] = fCells[iD].size(); + first = false; + } + fCells[iD].push_back(AliITSUCACell(fDoublets[iD][iD0].x,fDoublets[iD][iD0].y, + fDoublets[iD + 1][iD1].y,iD0,iD1,curv,n)); +#ifdef _TUNING_ + if (fGood) { + fTan->Fill(TMath::Abs(fDoublets[iD][iD0].tanL - fDoublets[iD + 1][iD1].tanL)); + fPhi->Fill(TMath::Abs(fDoublets[iD][iD0].phi - fDoublets[iD + 1][iD1].phi)); + fGDCAZ[iD]->Fill(fabs(extz-GetZ())); + numberOfGoodCells++; + } else { + fTanF->Fill(TMath::Abs(fDoublets[iD][iD0].tanL - fDoublets[iD + 1][iD1].tanL)); + fPhiF->Fill(TMath::Abs(fDoublets[iD][iD0].phi - fDoublets[iD + 1][iD1].phi)); + fFDCAZ[iD]->Fill(fabs(extz - GetZ())); + } + totalNumberOfCells++; +#endif + } + } + } + } + } + } + // + for (int iD = 0; iD < 4; ++iD) { + for (size_t c0 = 0; c0 < fCells[iD].size(); ++c0) { + const int idx = fCells[iD][c0].d1(); + for (size_t c1 = tLUT[iD][idx]; idx == fCells[iD + 1][c1].d0(); ++c1) { +#ifdef _TUNING_ + fGood = (fLayer[iD].GetClusterSorted(fCells[iD][c0].x())->GetLabel(0) == fLayer[iD + 1].GetClusterSorted(fCells[iD][c0].y())->GetLabel(0) && + fLayer[iD + 1].GetClusterSorted(fCells[iD][c0].y())->GetLabel(0) == fLayer[iD + 2].GetClusterSorted(fCells[iD][c0].z())->GetLabel(0) && + fLayer[iD + 2].GetClusterSorted(fCells[iD][c0].z())->GetLabel(0) == fLayer[iD + 3].GetClusterSorted(fCells[iD + 1][c1].z())->GetLabel(0) && + fLayer[iD].GetClusterSorted(fCells[iD][c0].x())->GetLabel(0) > 0); +#endif + float *n0 = fCells[iD][c0].GetN(); + float *n1 = fCells[iD + 1][c1].GetN(); + const float dn2 = ((n0[0] - n1[0]) * (n0[0] - n1[0]) + (n0[1] - n1[1]) * (n0[1] - n1[1]) + + (n0[2] - n1[2]) * (n0[2] - n1[2])); + const float dp = fabs(fCells[iD][c0].GetCurvature() - fCells[iD + 1][c1].GetCurvature()); + if (dn2 < kmaxDN[iD] && dp < kmaxDP[iD]) { +#ifdef _TUNING_ + assert(fCells[iD + 1][c1].Combine(fCells[iD][c0], c0)); + if (fGood) { + fGoodCombChi2[iD]->Fill(dp); + fGoodCombN[iD]->Fill(dn2); + cellsCombiningSuccesses++; + } else { + fFakeCombChi2[iD]->Fill(dp); + fFakeCombN[iD]->Fill(dn2); + cellsWrongCombinations++; + } +#else + fCells[iD + 1][c1].Combine(fCells[iD][c0], c0); +#endif + } + } + } + } + // + for (int level = 5; level > 1; --level) { + vector roads; + + for (int iCL = 4; iCL >= level - 1; --iCL) { + for (size_t iCell = 0; iCell < fCells[iCL].size(); ++iCell) { + if (fCells[iCL][iCell].GetLevel() != level) + { + continue; + } + roads.push_back(AliITSUCARoad(iCL,iCell)); + for(size_t iN = 0; iN < fCells[iCL][iCell].NumberOfNeighbours(); ++iN) { + const int currD = iCL - 1; + const int neigh = fCells[iCL][iCell](iN); + if(iN > 0) + { + roads.push_back(AliITSUCARoad(iCL,iCell)); + } + CellsTreeTraversal(roads,neigh,currD); + } + fCells[iCL][iCell].SetLevel(0u); // Level = -1 + } + } + + for (size_t iR = 0; iR < roads.size(); ++iR) + { + if (roads[iR].N != level) + continue; +#ifdef _TUNING_ + ++totalNumberOfRoads; +#endif + AliITSUTrackCooked tr; + int first = -1,last = -1; + ClsInfo_t *cl0 = 0x0,*cl1 = 0x0,*cl2 = 0x0; + for(int i = 0; i < 5; ++i) + { + if (roads[iR][i] < 0) + continue; + + if (first < 0) + { + cl0 = fLayer[i][fCells[i][roads[iR][i]].x()]; + tr.SetClusterIndex(i,fLayer[i][fCells[i][roads[iR][i]].x()]->index); + tr.SetClusterIndex(i + 1,fLayer[i + 1][fCells[i][roads[iR][i]].y()]->index); + first = i; + } + tr.SetClusterIndex(i + 2,fLayer[i + 2][fCells[i][roads[iR][i]].z()]->index); + last = i; + } + AliITSUClusterPix* c = fLayer[last + 2].GetClusterSorted(fCells[last][roads[iR][last]].z()); + cl2 = fLayer[last + 2][fCells[last][roads[iR][last]].z()]; + first = (last + first) / 2; +// if (roads[iR][first] < 0) { +// printf("roads[%lu][%i] : %i %i %i %i %i\n",iR,first,roads[iR][0],roads[iR][1],roads[iR][2],roads[iR][3],roads[iR][4]); +// } + cl1 = fLayer[first + 1][fCells[first][roads[iR][first]].y()]; + // Init track parameters + double cv = Curvature(cl0->x,cl0->y,cl1->x,cl1->y,cl2->x,cl2->y); + double tgl = TanLambda(cl0->x,cl0->y,cl2->x,cl2->y,cl0->z,cl2->z); + + AliITSUCATrackingStation::ITSDetInfo_t det = fLayer[last + 2].GetDetInfo(cl2->detid); + double x = det.xTF + c->GetX(); // I'd like to avoit using AliITSUClusterPix... + double alp = det.phiTF; + double par[5] = {c->GetY(),c->GetZ(),0,tgl,cv}; + double cov[15] = { + 5*5, + 0, 5*5, + 0, 0 , 0.7*0.7, + 0, 0, 0, 0.7*0.7, + 0, 0, 0, 0, 10 + }; + tr.Set(x,alp,par,cov); + AliITSUTrackCooked tt(tr); + tt.ResetClusters(); + if (RefitAt(2.1, &tt, &tr)) + new((*fCandidates[level - 2])[fCandidates[level - 2]->GetEntriesFast()]) AliITSUTrackCooked(tt); + } + } +#ifdef _TUNING_ +// Info("FindTracksCA","Iteration %i",iteration); + Info("FindTracksCA","Good doublets: %d",numberOfGoodDoublets); + Info("FindTracksCA","Number of doublets: %d",totalNumberOfDoublets); + Info("FindTracksCA","Good cells: %d",numberOfGoodCells); + Info("FindTracksCA","Number of cells: %d",totalNumberOfCells); +// Info("FindTracksCA","Cells combining inefficiencies: %d",cellsCombiningInefficiencies); + Info("FindTracksCA","Cells combining successes: %d",cellsCombiningSuccesses); + Info("FindTracksCA","Cells wrong combinations: %d",cellsWrongCombinations); +// Info("FindTracksCA","Number of found roads: %d",numberOfRoads); +// Info("FindTracksCA","Number of FULL found roads: %d",numberOfFullRoads); + Info("FindTracksCA","Roads survived after first selection: %d",totalNumberOfRoads); +#endif +} + +//__________________________________________________________________________________________________ +Int_t AliITSUCATracker::PropagateBack(AliESDEvent * event) +{ + // Here, we implement the Kalman smoother ? + // The clusters must be already loaded +// Int_t n=event->GetNumberOfTracks(); +// Int_t ntrk=0; +// Int_t ngood=0; +// for (Int_t i = 0; i < n; i++) +// { +// AliESDtrack *esdTrack=event->GetTrack(i); +// +// if ((esdTrack->GetStatus()&AliESDtrack::kITSin) == 0) +// continue; +// +// AliITSUTrackCooked track(*esdTrack); +// +// track.ResetCovariance(10.); +// +// int points[2 * AliITSUAux::kMaxLayers]; +// for (UInt_t k = 0; k < 2 * AliITSUAux::kMaxLayers; k++) +// points[k] = -1; +// Int_t nc = track.GetNumberOfClusters(); +// for (Int_t k = 0; k < nc; k++) +// { +// const int layer = (track.GetClusterIndex(k) & 0xf0000000) >> 28; +// const int idx = (track.GetClusterIndex(k) & 0x0fffffff); +// points[layer << 1]=idx; +// } +// +// if (RefitTrack(&track,points,40,1) >= 0) { +// +// CookLabel(&track, 0.); //For comparison only +// Int_t label = track.GetLabel(); +// if (label > 0) +// ngood++; +// +// esdTrack->UpdateTrackParams(&track,AliESDtrack::kITSout); +// ntrk++; +// } +// } +// +// Info("PropagateBack","Back propagated tracks: %d",ntrk); +// if (ntrk) +// Info("PropagateBack","Good tracks/back propagated: %f",Float_t(ngood)/ntrk); +// +// return 0; + + Int_t n=event->GetNumberOfTracks(); + Int_t ntrk=0; + Int_t ngood=0; + for (Int_t i=0; iGetTrack(i); + + if ((esdTrack->GetStatus()&AliESDtrack::kITSin)==0) continue; + + AliITSUTrackCooked track(*esdTrack); + AliITSUTrackCooked temp(*esdTrack); + + temp.ResetCovariance(10.); + temp.ResetClusters(); + + if (RefitAt(40., &temp, &track)) { + + CookLabel(&temp, 0.); //For comparison only + Int_t label = temp.GetLabel(); + if (label > 0) ngood++; + + esdTrack->UpdateTrackParams(&temp,AliESDtrack::kITSout); + //UseClusters(fTrackToFollow); + ntrk++; + } + } + + Info("PropagateBack","Back propagated tracks: %d",ntrk); + if (ntrk) + Info("PropagateBack","Good tracks/back propagated: %f",Float_t(ngood)/ntrk); + + return 0; +} + +//__________________________________________________________________________________________________ +Int_t AliITSUCATracker::RefitInward(AliESDEvent * event) +{ + // Some final refit, after the outliers get removed by the smoother ? + // The clusters must be loaded +// +// Int_t n = event->GetNumberOfTracks(); +// Int_t ntrk = 0; +// Int_t ngood = 0; +// for (Int_t i = 0; i < n; i++) { +// AliESDtrack *esdTrack=event->GetTrack(i); +// +// if ((esdTrack->GetStatus()&AliESDtrack::kITSout) == 0) continue; +// +// AliITSUTrackCooked track(*esdTrack); +// +// track.ResetCovariance(10.); +// +// int points[2 * AliITSUAux::kMaxLayers]; +// for (UInt_t k = 0; k < 2 * AliITSUAux::kMaxLayers; k++) +// points[k] = -1; +// Int_t nc = track.GetNumberOfClusters(); +// for (Int_t k = 0; k < nc; k++) +// { +// const int layer = (track.GetClusterIndex(k) & 0xf0000000) >> 28; +// const int idx = (track.GetClusterIndex(k) & 0x0fffffff); +// points[layer << 1] = idx; +// } +// +// if (RefitTrack(&track,points,1.8,1) >= 0) { //2.1,1)>=0) { +// +// //if (!track.PropagateTo(1.8, 2.27e-3, 35.28*1.848)) continue; +// CookLabel(&track, 0.); //For comparison only +// Int_t label=track.GetLabel(); +// if (label>0) ngood++; +// +// //cout << esdTrack->GetStatus() << " "; +// esdTrack->UpdateTrackParams(&track,AliESDtrack::kITSrefit); +// //cout << esdTrack->GetStatus() << endl; +// ntrk++; +// } +// } +// +// Info("RefitInward","Refitted tracks: %d",ntrk); +// if (ntrk) +// Info("RefitInward","Good tracks/refitted: %f",Float_t(ngood) / ntrk); +// +// return 0; + Int_t n = event->GetNumberOfTracks(); + Int_t ntrk = 0; + Int_t ngood = 0; + for (Int_t i = 0; i < n; i++) { + AliESDtrack *esdTrack = event->GetTrack(i); + + if ((esdTrack->GetStatus() & AliESDtrack::kITSout) == 0) continue; + + AliITSUTrackCooked track(*esdTrack); + AliITSUTrackCooked temp(*esdTrack); + + temp.ResetCovariance(10.); + temp.ResetClusters(); + + if (!RefitAt(2.1, &temp, &track)) continue; + //Cross the beam pipe + if (!temp.PropagateTo(1.8, 2.27e-3, 35.28 * 1.848)) continue; + + CookLabel(&temp, 0.); //For comparison only + Int_t label = temp.GetLabel(); + if (label > 0) ngood++; + + esdTrack->UpdateTrackParams(&temp,AliESDtrack::kITSrefit); + //esdTrack->RelateToVertex(event->GetVertex(),GetBz(),33.); + //UseClusters(fTrackToFollow); + ntrk++; + } + + Info("RefitInward","Refitted tracks: %d",ntrk); + if (ntrk) + Info("RefitInward","Good tracks/refitted: %f",Float_t(ngood)/ntrk); + + return 0; +} + +//__________________________________________________________________________________________________ +Int_t AliITSUCATracker::LoadClusters(TTree *cluTree) +{ + // This function reads the ITSU clusters from the tree, + // sort them, distribute over the internal tracker arrays, etc + + fITS->LoadClusters(cluTree); + fITS->ProcessClusters(); + // + // I consider a single vertex event for the moment. + double vertex[3] = {GetX(),GetY(),GetZ()}; + for (int iL = 0; iL < 7; ++iL) { + fLayer[iL].Init(fITS->GetLayerActive(iL),fITS->GetGeom()); + AliVertex v(vertex,1,1); + fLayer[iL].SortClusters(&v); + fUsedClusters[iL].resize(fLayer[iL].GetNClusters(),false); + } + return 0; +} + +//__________________________________________________________________________________________________ +void AliITSUCATracker::UnloadClusters() +{ + // This function unloads ITSU clusters from the memory + for (int i = 0;i < 7;++i) { + fLayer[i].Clear(); + fUsedClusters[i].clear(); + } + for (int i = 0; i < 6; ++i) { + fDoublets[i].clear(); + } + for (int i = 0; i < 5; ++i) { + fCells[i].clear(); + } + for (int i = 0; i < 4; ++i) + { + fCandidates[i]->Clear("C"); + } + // for(int i=0;i<5;++i) + // fCells[i].clear(); + +} + +//__________________________________________________________________________________________________ +Double_t AliITSUCATracker::RefitTrack(AliExternalTrackParam* trc, + Int_t clInfo[2*AliITSUAux::kMaxLayers], + Double_t rDest, Int_t stopCond) +{ + // refit track till radius rDest. + // if stopCond<0 : propagate till last cluster then stop + // if stopCond==0: propagate till last cluster then try to go till limiting + // rDest, don't mind if fail + // if stopCond>0 : rDest must be reached + // + // The clList should provide the indices of clusters at corresponding layer + // (as stored in the layer TClonesArray, with convention (allowing for up to 2 + // clusters per layer due to the overlaps): if there is a cluster on given + // layer I, then it should be stored at clInfo[2*I-1] if there is an + // additional cluster on this layer, it goes to clInfo[2*I], + // -1 means no cluster + // + double rCurr = Sqrt(trc->GetX()*trc->GetX() + trc->GetY()*trc->GetY()); + int dir,lrStart,lrStop; + // + dir = rCurrFindFirstLayerID(rCurr,dir); + lrStop = fITS->FindLastLayerID(rDest,dir);//lrid before which we have to stop + // + if (lrStop<0 || lrStart<0) + { + AliFatal(Form("Failed to find start(%d) or last(%d) layers. " + "Track from %.3f to %.3f",lrStart,lrStop,rCurr,rDest)); + } + // + int nCl = 0; + for (int i=2*fITS->GetNLayersActive();i--;) { + if (clInfo[i]<0) + continue; + nCl++; + } + // + AliExternalTrackParam tmpTr(*trc); + double chi2 = 0; + int iclLr[2],nclLr; + int nclFit = 0; + // + + int lrStop1 = lrStop+dir; + for (int ilr=lrStart;ilr!=lrStop1;ilr+=dir) { + AliITSURecoLayer* lr = fITS->GetLayer(ilr); + if ( dir*(rCurr-lr->GetR(dir))>0) continue; // this layer is already passed + int ilrA2,ilrA = lr->GetActiveID(); + // passive layer or active w/o hits will be traversed on the way to next + // cluster + if (!lr->IsActive() || clInfo[ilrA2=(ilrA<<1)]<0) continue; + // + // select the order in which possible 2 clusters (in case of the overlap) + // will be traversed and fitted + nclLr=0; + if (dir>0) { // clusters are stored in increasing radius order + iclLr[nclLr++]=clInfo[ilrA2++]; + if (clInfo[ilrA2]>=0) iclLr[nclLr++]=clInfo[ilrA2]; + } + else { + if ( clInfo[ilrA2+1]>=0 ) iclLr[nclLr++]=clInfo[ilrA2+1]; + iclLr[nclLr++]=clInfo[ilrA2]; + } + // + Bool_t transportedToLayer = kFALSE; + for (int icl=0;iclGetSensorFromID(clus->GetVolumeId()); + if (!tmpTr.Rotate(sens->GetPhiTF())) return -1; + // + double xClus = sens->GetXTF()+clus->GetX(); + if (!transportedToLayer) { + if (ilr!=lrStart && !TransportToLayerX(&tmpTr,lrStart,ilr,xClus)) { + return -1; // go to the entrance to the layer + } + lrStart = ilr; + transportedToLayer = kTRUE; + } + // + if (!PropagateSeed(&tmpTr,xClus,fCurrMass)) return -1; + // + Double_t p[2]={clus->GetY(), clus->GetZ()}; + Double_t cov[3]= \ + {clus->GetSigmaY2(), clus->GetSigmaYZ(), clus->GetSigmaZ2()}; + double chi2cl = tmpTr.GetPredictedChi2(p,cov); + chi2 += chi2cl; + // + if ( !tmpTr.Update(p,cov) ) return -1; + if (++nclFit==nCl && stopCond<0) { + *trc = tmpTr; + return chi2; // it was requested to not propagate after last update + } + } + // + } + // All clusters were succesfully fitted. Even if the track does not reach + // rDest, this is enough to validate it. Still, try to go as close as possible + // to rDest. + // + if (lrStart!=lrStop) { + if (!TransportToLayer(&tmpTr,lrStart,lrStop)) + return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory + if (!GoToExitFromLayer(&tmpTr,fITS->GetLayer(lrStop),dir)) + return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory + } + // go to the destination radius. Note that here we don't select direction + // to avoid precision problems + if (!tmpTr.GetXatLabR(rDest,rDest,GetBz(),0) || \ + !PropagateSeed(&tmpTr,rDest,fCurrMass, 100, kFALSE)) { + return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory + } + *trc = tmpTr; + + return chi2; +} + +//__________________________________________________________________________________________________ +Bool_t AliITSUCATracker::PropagateSeed(AliExternalTrackParam *seed, Double_t xToGo, Double_t mass, + Double_t maxStep, Bool_t matCorr) +{ + // propagate seed to given x applying material correction if requested + const Double_t kEpsilon = 1e-5; + Double_t xpos = seed->GetX(); + Int_t dir = (xpos < xToGo) ? 1 : -1; + Double_t xyz0[3],xyz1[3]; + // + Bool_t updTime = dir > 0 && seed->IsStartedTimeIntegral(); + if (matCorr || updTime) seed->GetXYZ(xyz1); //starting global position + while ( (xToGo - xpos) * dir > kEpsilon){ + Double_t step = dir * TMath::Min(TMath::Abs(xToGo-xpos), maxStep); + Double_t x = xpos + step; + Double_t bz = GetBz(); // getting the local Bz + if (!seed->PropagateTo(x,bz)) return kFALSE; + double ds = 0; + if (matCorr || updTime) { + xyz0[0] = xyz1[0]; // global pos at the beginning of step + xyz0[1] = xyz1[1]; + xyz0[2] = xyz1[2]; + seed->GetXYZ(xyz1); // // global pos at the end of step + // + if (matCorr) + { + Double_t xrho,xx0; + ds = GetMaterialBudget(xyz0,xyz1,xx0,xrho); + if (dir>0) xrho = -xrho; // outward should be negative + if (!seed->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE; + } + else + { // matCorr is not requested but time integral is + double d0 = xyz1[0] - xyz0[0]; + double d1 = xyz1[1] - xyz0[1]; + double d2 = xyz1[2] - xyz0[2]; + ds = TMath::Sqrt(d0 * d0 + d1 * d1 + d2 * d2); + } + } + if (updTime) seed->AddTimeStep(ds); + // + xpos = seed->GetX(); + } + return kTRUE; +} + +//__________________________________________________________________________________________________ +Bool_t AliITSUCATracker::TransportToLayer(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, + Double_t rLim) +{ + // transport track from layerFrom to the entrance of layerTo or to rLim (if>0), wathever is closer + // + if (lTo == lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo)); + // + int dir = lTo > lFrom ? 1 : -1; + // lrFr can be 0 when extrapolation from TPC to ITS is requested + AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); + Bool_t checkFirst = kTRUE; + Bool_t limReached = kFALSE; + while(lFrom != lTo) + { + if (lrFr) + { + if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) + { + return kFALSE; // go till the end of current layer + } + checkFirst = kFALSE; + } + AliITSURecoLayer* lrTo = fITS->GetLayer( (lFrom += dir) ); + if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom)); + // + // go the entrance of the layer, assuming no materials in between + double xToGo = lrTo->GetR(-dir); + if (rLim > 0) + { + if (dir > 0) + { + if (rLim < xToGo) + { + xToGo = rLim; + limReached = kTRUE; + } + } + else + { + if (rLim > xToGo) + { + xToGo = rLim; + limReached = kTRUE; + } + } + } + // double xts = xToGo; + if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) + { + // printf("FailHere1: %f %f %d\n",xts,xToGo,dir); + // seed->Print("etp"); + return kFALSE; + } + if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) + { + //printf("FailHere2: %f %f %d\n",xts,xToGo,dir); + //seed->Print("etp"); + return kFALSE; + } + lrFr = lrTo; + if (limReached) break; + } + return kTRUE; + // +} + +//__________________________________________________________________________________________________ +Bool_t AliITSUCATracker::TransportToLayerX(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, + Double_t xStop) +{ + // transport track from layerFrom to the entrance of layerTo but do not pass + // control parameter X + if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo)); + // + int dir = lTo > lFrom ? 1 : -1; + // lrFr can be 0 when extrapolation from TPC to ITS is requested + AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); + Bool_t checkFirst = kTRUE; + while(lFrom != lTo) + { + if (lrFr) { + if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) + { + return kFALSE; // go till the end of current layer + } + checkFirst = kFALSE; + } + AliITSURecoLayer* lrTo = fITS->GetLayer( (lFrom += dir) ); + if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom)); + // + // go the entrance of the layer, assuming no materials in between + double xToGo = lrTo->GetR(-dir); // R of the entrance to layer + // + // double xts = xToGo; + if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) { + // printf("FailHere1: %f %f %d\n",xts,xToGo,dir); + // seed->Print("etp"); + return kFALSE; + } + if ( (dir > 0 && xToGo > xStop) || (dir < 0 && xToGo < xStop) ) xToGo = xStop; + // + #ifdef _ITSU_DEBUG_ + AliDebug(2,Form("go in dir=%d to R=%.4f(X:%.4f)",dir,lrTo->GetR(-dir), xToGo)); + #endif + if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) { + //printf("FailHere2: %f %f %d\n",xts,xToGo,dir); + //seed->Print("etp"); + return kFALSE; + } + lrFr = lrTo; + } + return kTRUE; + // +} + +//__________________________________________________________________________________________________ +Bool_t AliITSUCATracker::GoToExitFromLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, + Int_t dir, Bool_t check) +{ + // go to the exit from lr in direction dir, applying material corrections in steps specific + // for this layer. + // If check is requested, do this only provided the track has not exited the layer already + double xToGo = lr->GetR(dir); + if (check) { // do we need to track till the surface of the current layer ? + double curR2 = seed->GetX() * seed->GetX() + seed->GetY() * seed->GetY(); // current radius + if (dir > 0) { + if (curR2 - xToGo*xToGo > -fgkToler) return kTRUE; + } // on the surface or outside of the layer + else if (dir < 0) { + if (xToGo * xToGo - curR2 > -fgkToler) return kTRUE; + } // on the surface or outside of the layer + } + if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE; + // go via layer to its boundary, applying material correction. + if (!PropagateSeed(seed,xToGo,fCurrMass, lr->GetMaxStep())) return kFALSE; + // + return kTRUE; + // +} + +//__________________________________________________________________________________________________ +Bool_t AliITSUCATracker::GoToEntranceToLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, + Int_t dir, Bool_t check) +{ + // go to the entrance of lr in direction dir, w/o applying material corrections. + // If check is requested, do this only provided the track did not reach the layer already + double xToGo = lr->GetR(-dir); + if (check) { // do we need to track till the surface of the current layer ? + double curR2 = seed->GetX() * seed->GetX() + seed->GetY() * seed->GetY(); // current radius + if (dir>0) { if (curR2 - xToGo * xToGo > -fgkToler) return kTRUE; } // already passed + else if (dir<0) { if (xToGo * xToGo - curR2 > -fgkToler) return kTRUE; } // already passed + } + if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE; + // go via layer to its boundary, applying material correction. + if (!PropagateSeed(seed,xToGo,fCurrMass, 100, kFALSE)) return kFALSE; + return kTRUE; + // +} + +//__________________________________________________________________________________________________ +Double_t AliITSUCATracker::GetMaterialBudget(const double* pnt0,const double* pnt1, double& x2x0, + double& rhol) const +{ + double par[7]; + if (fUseMatLUT && fMatLUT) { + double d = fMatLUT->GetMatBudget(pnt0,pnt1,par); + x2x0 = par[AliITSUMatLUT::kParX2X0]; + rhol = par[AliITSUMatLUT::kParRhoL]; + return d; + } + else { + MeanMaterialBudget(pnt0,pnt1,par); + x2x0 = par[1]; + rhol = par[0]*par[4]; + return par[4]; + } +} + +//__________________________________________________________________________________________________ +inline AliCluster* AliITSUCATracker::GetCluster(Int_t index) const +{ + const Int_t l=(index & 0xf0000000) >> 28; + const Int_t c=(index & 0x0fffffff); + return (AliCluster*)fLayer[l].GetClusterUnSorted(c) ; +} + +//__________________________________________________________________________________________________ +bool AliITSUCATracker::CellParams(int l, ClsInfo_t* c1, ClsInfo_t* c2, ClsInfo_t* c3, + float &curv, float n[3]) +{ + // Calculation of cell params and filtering using a DCA cut wrt beam line position. + + // Mapping the hits + const float mHit0[3] = {c1->x, c1->y, c1->r * c1->r}; + const float mHit1[3] = {c2->x, c2->y, c2->r * c2->r}; + const float mHit2[3] = {c3->x, c3->y, c3->r * c3->r}; + // Computing the deltas + const float mD10[3] = {mHit1[0] - mHit0[0],mHit1[1] - mHit0[1],mHit1[2] - mHit0[2]}; + const float mD20[3] = {mHit2[0] - mHit0[0],mHit2[1] - mHit0[1],mHit2[2] - mHit0[2]}; + // External product of the deltas -> n + n[0] = (mD10[1] * mD20[2]) - (mD10[2] * mD20[1]); + n[1] = (mD10[2] * mD20[0]) - (mD10[0] * mD20[2]); + n[2] = (mD10[0] * mD20[1]) - (mD10[1] * mD20[0]); + // Normalisation + const float norm = TMath::Sqrt((n[0] * n[0]) + (n[1] * n[1]) + (n[2] * n[2])); + if (norm < 1e-20f || fabs(n[2]) < 1e-20f) + return false; + n[0] /= norm; + n[1] /= norm; + n[2] /= norm; + // Center of the circle + const float c[2] = {-0.5f * n[0] / n[2], -0.5f * n[1] / n[2]}; + // Constant + const float k = - n[0] * mHit1[0] - n[1] * mHit1[1] - n[2] * mHit1[2]; + // Radius of the circle + curv = TMath::Sqrt((1.f - n[2] * n[2] - 4.f * k * n[2]) / (4.f * n[2] * n[2])); + // Distance of closest approach to the beam line + const float dca = fabs(curv - sqrt(c[0] * c[0] + c[1] * c[1])); + // Cut on the DCA + if (dca > kmaxDCAxy[l]) { + return false; + } +#ifdef _TUNING_ + if (fGood) { + fGDCAXY[l]->Fill(dca); + } else { + fFDCAXY[l]->Fill(dca); + } +#endif + + curv = 1.f / curv; + return true; +} + +//__________________________________________________________________________________________________ +Bool_t AliITSUCATracker::RefitAt(Double_t xx, AliITSUTrackCooked *t, const AliITSUTrackCooked *c) { + //-------------------------------------------------------------------- + // This function refits the track "t" at the position "x" using + // the clusters from "c" + //-------------------------------------------------------------------- + + const int nLayers = 7; + Int_t index[nLayers]; + Int_t k; + for (k = 0; k < nLayers; k++) index[k] = -1; + Int_t nc = c->GetNumberOfClusters(); + for (k = 0; k < nc; k++) { + Int_t idx = c->GetClusterIndex(k), nl = (idx&0xf0000000)>>28; + index[nl] = idx; + } + + Int_t from, to, step; + if (xx > t->GetX()) { + from = 0; + to = nLayers; + step = +1; + } else { + from = nLayers - 1; + to = -1; + step = -1; + } + + for (Int_t i = from; i != to; i += step) { + Int_t idx = index[i]; + if (idx >= 0) { + const AliCluster *cl = GetCluster(idx); + Float_t xr,ar; + cl->GetXAlphaRefPlane(xr, ar); + if (!t->Propagate(Double_t(ar), Double_t(xr), GetBz())) { + //Warning("RefitAt","propagation failed !\n"); + return kFALSE; + } + Double_t chi2 = t->GetPredictedChi2(cl); +// if (chi2 < 100) + t->Update(cl, chi2, idx); + } else { + Double_t r = fgkR[i]; + Double_t phi,z; + if (!t->GetPhiZat(r,phi,z)) { + //Warning("RefitAt","failed to estimate track !\n"); + return kFALSE; + } + if (!t->Propagate(phi, r, GetBz())) { + //Warning("RefitAt","propagation failed !\n"); + return kFALSE; + } + } + Double_t xx0 = (i > 2) ? 0.008 : 0.003; // Rough layer thickness + Double_t x0 = 9.36; // Radiation length of Si [cm] + Double_t rho = 2.33; // Density of Si [g/cm^3] + Double_t mass = t->GetMass(); + t->CorrectForMeanMaterial(xx0, - step * xx0 * x0 * rho, mass, kTRUE); + } + + if (!t->PropagateTo(xx,0.,0.)) return kFALSE; + return kTRUE; +} + diff --git a/ITS/UPGRADE/AliITSUCATracker.h b/ITS/UPGRADE/AliITSUCATracker.h new file mode 100644 index 00000000000..3c03d73dea6 --- /dev/null +++ b/ITS/UPGRADE/AliITSUCATracker.h @@ -0,0 +1,144 @@ +//-*- Mode: C++ -*- +// ************************************************************************ +// This file is property of and copyright by the ALICE ITSU Project * +// ALICE Experiment at CERN, All rights reserved. * +// See cxx source for full Copyright notice * +// * +//************************************************************************* + +#ifndef ALIITSUCATRACKER_H +#define ALIITSUCATRACKER_H + +#define _TUNING_ + +#include "AliTracker.h" +#include "AliITSUGeomTGeo.h" +#include +#include +#include "AliITSUMatLUT.h" +#include "AliITSUAux.h" +#include "AliExternalTrackParam.h" +#include "AliITSUCATrackingStation.h" +#include "AliITSUCACell.h" +#include "AliITSUTrackCooked.h" + +class AliITSUReconstructor; +class AliITSURecoDet; +class AliITSURecoLayer; +class AliITSUTrackCooked; + +class TTree; +class AliCluster; +class AliESDEvent; + +typedef struct AliITSUCATrackingStation::ClsInfo ClsInfo_t; + +using std::vector; + +//__________________________________________________________________________________________________ +class Doublets { +public: + Doublets(int xx = 0, int yy = 0, float tL = 0.f, float ph = 0.f) + : x((unsigned short)xx) + , y((unsigned short)yy) + , tanL(tL) + , phi(ph) {} + unsigned short x,y; + float tanL, phi; +}; + +//__________________________________________________________________________________________________ +class AliITSUCATracker : public AliTracker { +public: + AliITSUCATracker(AliITSUReconstructor* rec=0); + virtual ~AliITSUCATracker(); + + // These functions must be implemented + Int_t Clusters2Tracks(AliESDEvent *event); + Int_t PropagateBack(AliESDEvent *event); + Int_t RefitInward(AliESDEvent *event); + Int_t LoadClusters(TTree *ct); + void UnloadClusters(); + AliCluster *GetCluster(Int_t index) const; + + // Possibly, other public functions + void Init(AliITSUReconstructor* rec); + Double_t RefitTrack(AliExternalTrackParam* trc, Int_t clInfo[2 * AliITSUAux::kMaxLayers], Double_t rDest, Int_t stopCond); + Bool_t PropagateSeed(AliExternalTrackParam *seed, Double_t xToGo, Double_t mass, Double_t maxStep=1.0, Bool_t matCorr=kTRUE); + Double_t GetMaterialBudget(const double* pnt0, const double* pnt1, double& x2x0, double& rhol) const; + Bool_t GoToEntranceToLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check=kFALSE); + Bool_t GoToExitFromLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check=kTRUE); + Bool_t TransportToLayerX(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t xStop); + Bool_t TransportToLayer(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t rLim=-1); + + void SetChi2Cut(float cut) { fChi2Cut = cut; } + void SetPhiCut(float cut) { fPhiCut = cut; } + void SetZCut(float cut) { fZCut = cut; } + +#ifdef _TUNING_ + bool fGood; + TH1F * fGoodCombChi2[5]; + TH1F * fFakeCombChi2[5]; + TH1F * fGoodCombN[4]; + TH1F * fFakeCombN[4]; + TH1F * fGDZ[6]; + TH1F * fGDXY[6]; + TH1F * fFDZ[6]; + TH1F * fFDXY[6]; + TH1F * fGDCAZ[5]; + TH1F * fGDCAXY[5]; + TH1F * fFDCAZ[5]; + TH1F * fFDCAXY[5]; + TH1F * fTan; + TH1F * fTanF; + TH1F * fPhi; + TH1F * fPhiF; + TH1F * fNEntries; +#endif + // +protected: + AliITSUCATracker(const AliITSUCATracker&); + + // void MakeTriplets(); + void CellsTreeTraversal(vector &roads, const int &iD, const int &doubl); + Bool_t InitTrackParams(AliITSUTrackCooked &track, int points[]); + void FindTracksCA(int iteration); + void GlobalFit(); + void MergeTracks( vector &vec, bool flags[] ); + float FilterSeed(AliITSUCACell &c1, AliITSUCACell &c2, int lrStart); + bool CellParams(int l, ClsInfo_t* c1, ClsInfo_t* c2, ClsInfo_t* c3, float &curv, float np[3]); + + Bool_t RefitAt(Double_t xx, AliITSUTrackCooked *t, const AliITSUTrackCooked *c); +private: + AliITSUCATracker &operator=(const AliITSUCATracker &tr); + void SetLabel(AliITSUTrackCooked &t, Float_t wrong); + + // Data members + + // classes for interfacing the geometry, materials etc. + AliITSUReconstructor* fReconstructor; // ITS global reconstructor + AliITSURecoDet* fITS; // interface to ITS, borrowed from reconstructor + AliITSUMatLUT* fMatLUT; // material lookup table + Bool_t fUseMatLUT; //! use material lookup table rather than TGeo + Double_t fCurrMass; // assumption about particle mass + + + // Internal tracker arrays, layers, modules, etc + AliITSUCATrackingStation fLayer[7]; + vector fUsedClusters[7]; + Float_t fChi2Cut; + Float_t fPhiCut; + Float_t fZCut; + vector fDoublets[6]; + vector fCells[5]; + TClonesArray *fCandidates[4]; + + static const Double_t fgkToler; // tracking tolerance + static const Double_t fgkChi2Cut; // chi2 cut during track merging + static const int fgkNumberOfIterations; + static const float fgkR[7]; + // + ClassDef(AliITSUCATracker,2) //ITSU stand-alone tracker +}; + +#endif // ALIITSUCATRACKER_H diff --git a/ITS/UPGRADE/AliITSUCATrackingStation.cxx b/ITS/UPGRADE/AliITSUCATrackingStation.cxx new file mode 100644 index 00000000000..a8ff3ef5686 --- /dev/null +++ b/ITS/UPGRADE/AliITSUCATrackingStation.cxx @@ -0,0 +1,368 @@ +//-*- Mode: C++ -*- +// ************************************************************************** +// This file is property of and copyright by the ALICE ITSU Project * +// ALICE Experiment at CERN, All rights reserved. * +// * +// Primary Author: Ruben Shahoyan * +// * +// Adapted to ITSU: Maximiliano Puccio * +// for the ITS Upgrade project * +// * +// Permission to use, copy, modify and distribute this software and its * +// documentation strictly for non-commercial purposes is hereby granted * +// without fee, provided that the above copyright notice appears in all * +// copies and that both the copyright notice and this permission notice * +// appear in the supporting documentation. The authors make no claims * +// about the suitability of this software for any purpose. It is * +// provided "as is" without express or implied warranty. * +// * +//*************************************************************************** + +#include "AliITSUCATrackingStation.h" +#include +#include "AliITSUGeomTGeo.h" +#include "AliVertex.h" +#include +#include +#include +#include "AliITSUAux.h" +#include "AliITSURecoSens.h" +#include + +using AliITSUAux::BringTo02Pi; + +//_________________________________________________________________ +AliITSUCATrackingStation::AliITSUCATrackingStation() +:fClusters(0x0) +,fVIDOffset(0) +,fNClusters(0) +,fZMin(0) +,fZMax(0) +,fDZInv(-1) +,fDPhiInv(-1) +,fNZBins(20) +,fNPhiBins(20) +,fQueryZBmin(-1) +,fQueryZBmax(-1) +,fQueryPhiBmin(-1) +,fQueryPhiBmax(-1) +,fBins(0) +,fOccBins(0) +,fNOccBins(0) +,fNFoundClusters(0) +,fFoundClusterIterator(0) +,fFoundBinIterator(0) +,fIndex() +,fFoundBins(0) +,fSortedClInfo(0) +{ + // def. c-tor +} + +//_________________________________________________________________ +AliITSUCATrackingStation::AliITSUCATrackingStation(int nzbins,int nphibins, + AliITSURecoLayer *lr, AliITSUGeomTGeo *geo) +:fClusters(0x0) +,fVIDOffset() +,fNClusters(0) +,fZMin(0.f) +,fZMax(0.f) +,fDZInv(-1) +,fDPhiInv(-1) +,fNZBins(nzbins) +,fNPhiBins(nphibins) +,fQueryZBmin(-1) +,fQueryZBmax(-1) +,fQueryPhiBmin(-1) +,fQueryPhiBmax(-1) +,fBins(0) +,fOccBins(0) +,fNOccBins(0) +,fNFoundClusters(0) +,fFoundClusterIterator(0) +,fFoundBinIterator(0) +,fIndex() +,fFoundBins(0) +,fSortedClInfo(0) +{ + // c-tor + Init(lr,geo); +} + +//_________________________________________________________________ +AliITSUCATrackingStation::~AliITSUCATrackingStation() +{ + // d-tor + delete[] fBins; + delete[] fOccBins; +// delete fClusters; +} + +//_________________________________________________________________ +void AliITSUCATrackingStation::Init(AliITSURecoLayer *lr, AliITSUGeomTGeo *geo) +{ + if (fClusters) { + printf("Already initialized\n"); + return; + } + + fClusters = *lr->GetClustersAddress(); + fZMin = lr->GetZMin(); + fZMax = lr->GetZMax(); + if (fNZBins < 1) fNZBins = 2; + if (fNPhiBins < 1) fNPhiBins = 1; + fDZInv = fNZBins / (fZMax - fZMin); + fDPhiInv = fNPhiBins / TMath::TwoPi(); + // + fBins = new ClBinInfo_t[fNZBins * fNPhiBins]; + fOccBins = new int[fNZBins * fNPhiBins]; + fNClusters = fClusters->GetEntriesFast(); + fSortedClInfo.reserve(fClusters->GetEntriesFast()); + fVIDOffset = ((AliITSUClusterPix*)fClusters->UncheckedAt(0))->GetVolumeId(); + // + // prepare detectors info + int detID = -1; + fIndex.resize(lr->GetNSensors(),-1); + fDetectors.reserve(lr->GetNSensors()); + for (int iCl = 0; iCl < fClusters->GetEntriesFast(); ++iCl) + { //Fill this layer with detectors + AliITSUClusterPix* c = (AliITSUClusterPix*)fClusters->UncheckedAt(iCl); + if (detID == c->GetVolumeId()) + { + continue; + } + detID = c->GetVolumeId(); + ITSDetInfo_t det; + det.index = iCl; + // + TGeoHMatrix m; + geo->GetOrigMatrix(detID,m); + // + fIndex[detID - fVIDOffset] = fDetectors.size(); + AliITSURecoSens* sens = lr->GetSensorFromID(detID); + const TGeoHMatrix *tm = geo->GetMatrixT2L(detID); + m.Multiply(tm); + double txyz[3] = {0.,0.,0.}, xyz[3] = {0.,0.,0.}; + m.LocalToMaster(txyz,xyz); + det.xTF = sens->GetXTF(); // TMath::Sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]); + + det.phiTF = sens->GetPhiTF();//TMath::ATan2(xyz[1],xyz[0]); + det.sinTF = TMath::Sin(det.phiTF); + det.cosTF = TMath::Cos(det.phiTF); + // + // compute the real radius (with misalignment) + TGeoHMatrix mmisal(*(geo->GetMatrix(detID))); + mmisal.Multiply(tm); + xyz[0] = 0.; + xyz[1] = 0.; + xyz[2] = 0.; + mmisal.LocalToMaster(txyz,xyz); + det.xTFmisal = TMath::Sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]); +// if (!TMath::AreEqualAbs(det.xTFmisal, sens->GetXTF(), 1e-9)) { +// cout << "Error: " << lr->GetActiveID() << " " << detID << " " << det.xTFmisal - sens->GetXTF()<< endl; +// } + + fDetectors.push_back(det); + } // end loop on detectors +} + +//_________________________________________________________________ +void AliITSUCATrackingStation::SortClusters(const AliVertex* vtx) +{ + // sort clusters and build fast lookup table + // + ClearSortedInfo(); + fSortedClInfo.reserve(fNClusters); + // + ClsInfo_t cl; + for (int icl = fNClusters;icl--;) { + AliITSUClusterPix* cluster = (AliITSUClusterPix*)fClusters->UncheckedAt(icl); + cluster->GetGlobalXYZ( (float*)&cl ); + // + if (vtx) { // phi and r will be computed wrt vertex + cl.x -= vtx->GetX(); + cl.y -= vtx->GetY(); + } + // + cl.r = TMath::Sqrt(cl.x*cl.x + cl.y*cl.y); + cl.phi = TMath::ATan2(cl.y,cl.x); + BringTo02Pi(cl.phi); + cl.index = icl; + cl.zphibin = GetBinIndex(GetZBin(cl.z),GetPhiBin(cl.phi)); + cl.detid = cluster->GetVolumeId() - fVIDOffset; + // + fSortedClInfo.push_back( cl ); + // + } + sort(fSortedClInfo.begin(), fSortedClInfo.end()); // sort in phi, z + // + // fill cells in phi,z + int currBin = -1; + for (int icl = 0;icl < fNClusters; ++icl) + { + ClsInfo_t &t = fSortedClInfo[icl]; + if (t.zphibin > currBin) + { // register new occupied bin + currBin = t.zphibin; + fBins[currBin].first = icl; + fBins[currBin].index = fNOccBins; + fOccBins[fNOccBins++] = currBin; + } + fBins[currBin].ncl++; + } + // Print("clb"); //RS +} + +//_________________________________________________________________ +void AliITSUCATrackingStation::Clear(Option_t *) +{ + // clear cluster info + ClearSortedInfo(); + fIndex.clear(); + fNClusters = 0; + if (fClusters) fClusters = 0x0; + // +} + +//_________________________________________________________________ +void AliITSUCATrackingStation::ClearSortedInfo() +{ + // clear cluster info + fSortedClInfo.clear(); + memset(fBins,0,fNZBins * fNPhiBins * sizeof(ClBinInfo_t)); + memset(fOccBins,0,fNZBins * fNPhiBins * sizeof(int)); + fNOccBins = 0; +} + +//_________________________________________________________________ +void AliITSUCATrackingStation::Print(Option_t *opt) const +{ + // dump cluster bins info + TString opts = opt; + opts.ToLower(); + printf("Stored %d clusters in %d occupied bins\n",fNClusters,fNOccBins); + // + if (opts.Contains("c")) { + printf("\nCluster info\n"); + for (int i = 0; i < fNClusters;i++) + { + const ClsInfo_t &t = fSortedClInfo[i]; + printf("#%5d Bin(phi/z):%03d/%03d Z:%+8.3f Phi:%+6.3f R:%7.3f Ind:%d ", + i,t.zphibin/fNZBins,t.zphibin%fNZBins,t.z,t.phi,t.r,t.index); + if (opts.Contains("l")) { // mc labels + AliITSUClusterPix* rp = (AliITSUClusterPix*)fClusters->UncheckedAt(t.index); + for (int l = 0;l < 3; l++) if (rp->GetLabel(l) >= 0) printf("| %d ",rp->GetLabel(l)); + } + printf("\n"); + } + } + // + if (opts.Contains("b")) { + printf("\nBins info (occupied only)\n"); + for (int i=0;i %3d cl from %d\n",i,fOccBins[i],fOccBins[i]/fNZBins,fOccBins[i]%fNZBins, + fBins[fOccBins[i]].ncl,fBins[fOccBins[i]].first); + } + } + // +} + +//_____________________________________________________________ +int AliITSUCATrackingStation::SelectClusters(float zmin,float zmax,float phimin,float phimax) +{ + // prepare occupied bins in the requested region + //printf("Select: Z %f %f | Phi: %f %f\n",zmin,zmax,phimin,phimax); + if (!fNOccBins) return 0; + if (zmax < fZMin || zmin > fZMax || zmin > zmax) return 0; + fFoundBins.clear(); + + fQueryZBmin = GetZBin(zmin); + if (fQueryZBmin < 0) fQueryZBmin = 0; + fQueryZBmax = GetZBin(zmax); + if (fQueryZBmax >= fNZBins) fQueryZBmax = fNZBins - 1; + BringTo02Pi(phimin); + BringTo02Pi(phimax); + fQueryPhiBmin = GetPhiBin(phimin); + fQueryPhiBmax = GetPhiBin(phimax); + int dbz = 0; + fNFoundClusters = 0; + int nbcheck = fQueryPhiBmax - fQueryPhiBmin + 1; //TODO:(MP) check if a circular buffer is feasible + if (nbcheck > 0) + { // no wrapping around 0-2pi, fast case + for (int ip = fQueryPhiBmin;ip <= fQueryPhiBmax;ip++) { + int binID = GetBinIndex(fQueryZBmin,ip); + if ( !(dbz = (fQueryZBmax-fQueryZBmin)) ) { // just one Z bin in the query range + ClBinInfo_t& binInfo = fBins[binID]; + if (!binInfo.ncl) continue; + fNFoundClusters += binInfo.ncl; + fFoundBins.push_back(binID); + continue; + } + int binMax = binID + dbz; + for ( ; binID <= binMax; binID++) { + ClBinInfo_t& binInfo = fBins[binID]; + if (!binInfo.ncl) continue; + fNFoundClusters += binInfo.ncl; + fFoundBins.push_back(binID); + } + } + } + else + { // wrapping + nbcheck += fNPhiBins; + for (int ip0 = 0;ip0 <= nbcheck;ip0++) { + int ip = fQueryPhiBmin + ip0; + if (ip >= fNPhiBins) ip -= fNPhiBins; + int binID = GetBinIndex(fQueryZBmin,ip); + if ( !(dbz = (fQueryZBmax - fQueryZBmin)) ) { // just one Z bin in the query range + ClBinInfo_t& binInfo = fBins[binID]; + if (!binInfo.ncl) continue; + fNFoundClusters += binInfo.ncl; + fFoundBins.push_back(binID); + continue; + } + int binMax = binID + dbz; + for (;binID <= binMax;binID++) { + ClBinInfo_t& binInfo = fBins[binID]; + if (!binInfo.ncl) continue; + fNFoundClusters += binInfo.ncl; + fFoundBins.push_back(binID); + } + } + } + fFoundClusterIterator = fFoundBinIterator = 0; + /* + //printf("Selected -> %d cl in %d bins\n",fNFoundClusters,(int)fFoundBins.size()); + for (int i=0;i<(int)fFoundBins.size();i++) { + int bn = fFoundBins[i]; + ClBinInfo_t& bin=fBins[bn]; + printf("#%d b:%d 1st: %3d Ncl:%d\n",i,bn,bin.first,bin.ncl); + } + printf("\n"); + */ + return fNFoundClusters; +} + +//_____________________________________________________________ +int AliITSUCATrackingStation::GetNextClusterInfoID() +{ + if (fFoundBinIterator < 0) return 0; + int currBin = fFoundBins[fFoundBinIterator]; + if (fFoundClusterIterator < fBins[currBin].ncl) { // same bin + return fBins[currBin].first + fFoundClusterIterator++; + } + if (++fFoundBinIterator < int(fFoundBins.size())) { // need to change bin + currBin = fFoundBins[fFoundBinIterator]; + fFoundClusterIterator = 1; + return fBins[currBin].first; + } + fFoundBinIterator = -1; + return -1; +} + +//_____________________________________________________________ +void AliITSUCATrackingStation::ResetFoundIterator() +{ + // prepare for a new loop over found clusters + if (fNFoundClusters) fFoundClusterIterator = fFoundBinIterator = 0; +} diff --git a/ITS/UPGRADE/AliITSUCATrackingStation.h b/ITS/UPGRADE/AliITSUCATrackingStation.h new file mode 100644 index 00000000000..3c54ac8919a --- /dev/null +++ b/ITS/UPGRADE/AliITSUCATrackingStation.h @@ -0,0 +1,159 @@ +//-*- Mode: C++ -*- +// ************************************************************************ +// This file is property of and copyright by the ALICE ITSU Project * +// ALICE Experiment at CERN, All rights reserved. * +// See cxx source for full Copyright notice * +// * +//************************************************************************* + +#ifndef ALIITSUCATRACKINGSTATION_H +#define ALIITSUCATRACKINGSTATION_H + +#include +#include +#include +#include +#include "AliITSURecoLayer.h" +#include "AliITSUClusterPix.h" +#include + +class AliVertex; + +class AliITSUCATrackingStation +{ + +public: + struct ClsInfo // cluster info, optionally XY origin at vertex + { + float x,y,z,phi,r; // lab params + int zphibin; // bins is z,phi + int index; // index in RecPoints array + int detid; // detector index //RS ??? Do we need it? + bool operator<(const ClsInfo &rhs) const {return zphibin < rhs.zphibin;} + // + }; + typedef struct ClsInfo ClsInfo_t; + // + struct ClBinInfo // info on bin clusters start, number of clusters + { + unsigned short ncl; // number of clusters + unsigned short first; // entry of 1st cluster in sorted vector of ClsInfo + int index; // index in the vector containing cells with non-0 occupancy + }; + typedef struct ClBinInfo ClBinInfo_t; + // + struct ITSDetInfo // info on sensor + { + int index; // sensor vid + float xTF,xTFmisal,phiTF,sinTF,cosTF; //tracking frame parameters of the detector + }; + typedef struct ITSDetInfo ITSDetInfo_t; + + + AliITSUCATrackingStation(); + AliITSUCATrackingStation(int nzbins,int nphibins, AliITSURecoLayer *lr, AliITSUGeomTGeo *geo); + virtual ~AliITSUCATrackingStation(); + AliITSUCATrackingStation::ClsInfo_t* operator[](int i) const + {return (AliITSUCATrackingStation::ClsInfo_t*)&fSortedClInfo[i];} + // + int GetVIDOffset() const {return fVIDOffset;} + int GetNClusters() const {return fNClusters;} + int GetNZBins() const {return fNZBins;} + int GetNPhiBins() const {return fNPhiBins;} + float GetZMin() const {return fZMin;} + float GetZMax() const {return fZMax;} + // + void SetNZBins(int v) {fNZBins = v;} + void SetNPhiBins(int v) {fNPhiBins = v;} + void SetZMin(float v) {fZMin = v;} + void SetZMax(float v) {fZMax = v;} + // + void Init(AliITSURecoLayer *lr, AliITSUGeomTGeo *geo); + // + void SortClusters(const AliVertex* vtx = 0); + int GetPhiBin(float phi) const {return phi * fDPhiInv;} + int GetZBin (float z) const {return (z - fZMin) * fDZInv;} + int GetBinIndex(int iz, int iphi) const {return iphi * fNZBins + iz;} + int GetBinZ(int ipz) const {return ipz % fNZBins;} + int GetBinPhi(int ipz) const {return ipz / fNZBins;} + void GetBinZPhi(int ipz,int &iz,int &iphi) const {iz = GetBinZ(ipz); iphi=GetBinPhi(ipz);} + // + int SelectClusters(float zmin,float zmax,float phimin,float phimax); + int GetNFoundBins() const {return fFoundBins.size();} + int GetFoundBin(int i) const {return fFoundBins[i];} + int GetFoundBinClusters(int i, int &first) const; + void ResetFoundIterator(); + AliITSUCATrackingStation::ClsInfo_t* GetClusterInfo(int i) const + {return (AliITSUCATrackingStation::ClsInfo_t*)&fSortedClInfo[i];} + AliITSUCATrackingStation::ClsInfo_t* GetNextClusterInfo(); + int GetNextClusterInfoID(); + AliITSUClusterPix* GetNextCluster(); + AliITSUClusterPix* GetClusterSorted(int i) const + {return (AliITSUClusterPix*)fClusters->UncheckedAt(fSortedClInfo[i].index);} + AliITSUClusterPix* GetClusterUnSorted(int i) const + {return (AliITSUClusterPix*)fClusters->UncheckedAt(i);} + // + AliITSUCATrackingStation::ITSDetInfo_t& GetDetInfo(int id) const + {assert(fIndex[id] > -1 && "Empty sensor");return (ITSDetInfo_t&)fDetectors[fIndex[id]];} + int GetNDetectors() const + {return fDetectors.size();} + + void ClearSortedInfo(); + virtual void Clear(Option_t *opt=""); + virtual void Print(Option_t *opt="") const; + +protected: + TClonesArray* fClusters; // externally supplied clusters + AliITSURecoLayer* fLayer; // layer id + int fVIDOffset; // offset of VID for detectors of this layer + int fNClusters; // N clusters + // + float fZMin; // Zmin + float fZMax; // Zmax + float fDZInv; // inverse size of Z bin + float fDPhiInv; // inverse size of Phi bin + int fNZBins; // N cells in Z + int fNPhiBins; // N cells in Phi + // + int fQueryZBmin; // min bin in Z of the query + int fQueryZBmax; // max bin in Z of the query + int fQueryPhiBmin; // min bin in phi of the query + int fQueryPhiBmax; // max bin in phi of the query + ClBinInfo_t* fBins; // 2D (z,phi) grid of clusters binned in z,phi + int* fOccBins; // id's of bins with non-0 occupancy + int fNOccBins; // number of occupied bins + int fNFoundClusters; // number of found clusters in the query zone + int fFoundClusterIterator; // at which cluster within the bin we are? + int fFoundBinIterator; // at which foune bin we are? + std::vector fIndex; + std::vector fFoundBins; // occupied bins satisfying to query + std::vector fSortedClInfo; // processed cluster info + std::vector fDetectors; // detector params + // +}; + +//_____________________________________________________ +inline int AliITSUCATrackingStation::GetFoundBinClusters(int i, int &first) const { + // set the entry of the first cl.info in the fSortedClInfo + // and return n clusters in the bin + ClBinInfo_t& bin = fBins[GetFoundBin(i)]; + first = bin.first; + return bin.ncl; +} + +//_____________________________________________________ +inline AliITSUClusterPix* AliITSUCATrackingStation::GetNextCluster() { + // return next cluster + ClsInfo_t* cli = GetNextClusterInfo(); + return cli ? (AliITSUClusterPix*)fClusters->UncheckedAt(cli->index) : 0; +} + +//_____________________________________________________________ +inline AliITSUCATrackingStation::ClsInfo_t* AliITSUCATrackingStation::GetNextClusterInfo() +{ + // return cluster info for next matching cluster + int id = GetNextClusterInfoID(); + return id < 0 ? 0 : (AliITSUCATrackingStation::ClsInfo_t*)&fSortedClInfo[id]; +} + +#endif diff --git a/ITS/UPGRADE/AliITSUTrackerSA.cxx b/ITS/UPGRADE/AliITSUTrackerSA.cxx deleted file mode 100644 index 9bd63b68be7..00000000000 --- a/ITS/UPGRADE/AliITSUTrackerSA.cxx +++ /dev/null @@ -1,885 +0,0 @@ -//-------------------------------------------------------------------------------- -// Implementation of the ITS tracker class -// It reads AliITSUClusterPix clusters and and fills the ESD with tracks -// -// The algorithm implemented here takes inspiration from UniCA code of FIAS -// group. -//-------------------------------------------------------------------------------- - -#include -#include -using TMath::Abs; -using TMath::Sort; -using TMath::Sqrt; -#include -#include -using std::sort; - -// Vc library -//#include "Vc/Vc" -//#include "AliITSUTrackerSAauxVc.h" // Structs and other stuff using Vc library -#include "AliLog.h" -#include "AliESDEvent.h" -#include "AliITSUClusterPix.h" -#include "AliITSUTrackerSA.h" -#include "AliITSUReconstructor.h" -#include "AliITSURecoDet.h" -#include "AliESDtrack.h" - -#include - -using std::cout; -using std::endl; -using std::flush; - -//#include "AliITSUtrackSA.h" // Some dedicated SA track class ? - -ClassImp(AliITSUTrackerSA) - -const Double_t AliITSUTrackerSA::fgkToler = 1e-6;// tolerance for layer on-surface check -const Double_t AliITSUTrackerSA::fgkChi2Cut = 600.f; - -//________________________________________________________________________________ -AliITSUTrackerSA::AliITSUTrackerSA(AliITSUReconstructor* rec) : -fReconstructor(rec), -fITS(0), -fMatLUT(0), -fUseMatLUT(kFALSE), -fCurrMass(0.14), -// -fClustersTC(), -fChi2Cut( fgkChi2Cut ), -fPhiCut( 1 ), -fRPhiCut( 1 ), -fZCut( 1 ) -{ - //-------------------------------------------------------------------- - // This default constructor needs to be provided - //-------------------------------------------------------------------- - if (rec) Init(rec); -} - -//________________________________________________________________________________ -AliITSUTrackerSA::AliITSUTrackerSA(const AliITSUTrackerSA &t): AliTracker(t), - fReconstructor(t.fReconstructor), - fITS(t.fITS), - fMatLUT(t.fMatLUT), - fUseMatLUT(t.fUseMatLUT), - fCurrMass(t.fCurrMass), - // - fClustersTC(), - fChi2Cut(fgkChi2Cut), - fPhiCut(), - fRPhiCut(), - fZCut() -{ - //-------------------------------------------------------------------- - // The copy constructor is protected - //-------------------------------------------------------------------- -} - -//________________________________________________________________________________ -AliITSUTrackerSA::~AliITSUTrackerSA() -{ - // d-tor - delete fMatLUT; -} - - -//_________________________________________________________________________ -void AliITSUTrackerSA::Init(AliITSUReconstructor* rec) -{ - // init with external reconstructor - // - fITS = rec->GetITSInterface(); - // - // create material lookup table - const int kNTest = 1000; - const double kStepsPerCM=5; - fMatLUT = new AliITSUMatLUT(fITS->GetRMin(),fITS->GetRMax(),Nint(kStepsPerCM*(fITS->GetRMax()-fITS->GetRMin()))); - double zmn = 1e6; - for (int ilr=fITS->GetNLayers();ilr--;) { - AliITSURecoLayer* lr = fITS->GetLayer(ilr); - if (zmn>Abs(lr->GetZMin())) zmn = Abs(lr->GetZMin()); - if (zmn>Abs(lr->GetZMax())) zmn = Abs(lr->GetZMax()); - } - fMatLUT->FillData(kNTest,-zmn,zmn); - // -} - -//________________________________________________________________________________ -Int_t AliITSUTrackerSA::Clusters2Tracks(AliESDEvent *event) { - //-------------------------------------------------------------------- - // This is the main tracking function - // The clusters must already be loaded - //-------------------------------------------------------------------- - - // Possibly, create the track "seeds" (combinatorial) - - // Possibly, increment the seeds with additional clusters (Kalman) - - // Possibly, (re)fit the found tracks - - // Three iterations: - // - High momentum first; - // - Low momentum with vertex constraint; - // - Everything else. - - CellsCreation(0); - CellularAutomaton(event); - // VertexFinding(); - // CellsCreation(1); - // CellularAutomaton(event); - - return 0; -} - -//________________________________________________________________________________ -Int_t AliITSUTrackerSA::PropagateBack(AliESDEvent * event) { - //-------------------------------------------------------------------- - // Here, we implement the Kalman smoother ? - // The clusters must be already loaded - //-------------------------------------------------------------------- - Int_t n=event->GetNumberOfTracks(); - Int_t ntrk=0; - Int_t ngood=0; - for (Int_t i=0; iGetTrack(i); - - if ((esdTrack->GetStatus()&AliESDtrack::kITSin)==0) continue; - - AliITSUTrackCooked track(*esdTrack); - - track.ResetCovariance(10.); - - int points[2*AliITSUAux::kMaxLayers]; - for (UInt_t k=0; k<2*AliITSUAux::kMaxLayers; k++) - points[k]=-1; - Int_t nc=track.GetNumberOfClusters(); - for (Int_t k=0; k>28; - const int idx = (track.GetClusterIndex(k)&0x0fffffff); - points[layer<<1]=idx; - } - - if (RefitTrack(&track,points,40,1)>=0) { - - CookLabel(&track, 0.); //For comparison only - Int_t label=track.GetLabel(); - if (label>0) ngood++; - - esdTrack->UpdateTrackParams(&track,AliESDtrack::kITSout); - ntrk++; - } - } - - Info("PropagateBack","Back propagated tracks: %d",ntrk); - if (ntrk) - Info("PropagateBack","Good tracks/back propagated: %f",Float_t(ngood)/ntrk); - - return 0; -} - -//________________________________________________________________________________ -Int_t AliITSUTrackerSA::RefitInward(AliESDEvent * event) { - //-------------------------------------------------------------------- - // Some final refit, after the outliers get removed by the smoother ? - // The clusters must be loaded - //-------------------------------------------------------------------- - Int_t n=event->GetNumberOfTracks(); - Int_t ntrk=0; - Int_t ngood=0; - for (Int_t i=0; iGetTrack(i); - - if ((esdTrack->GetStatus()&AliESDtrack::kITSout)==0) continue; - - AliITSUTrackCooked track(*esdTrack); - - track.ResetCovariance(10.); - - int points[2*AliITSUAux::kMaxLayers]; - for (UInt_t k=0; k<2*AliITSUAux::kMaxLayers; k++) - points[k]=-1; - Int_t nc=track.GetNumberOfClusters(); - for (Int_t k=0; k>28; - const int idx = (track.GetClusterIndex(k)&0x0fffffff); - points[layer<<1]=idx; - } - - if (RefitTrack(&track,points,1.8,1)>=0) { //2.1,1)>=0) { - - //if (!track.PropagateTo(1.8, 2.27e-3, 35.28*1.848)) continue; - CookLabel(&track, 0.); //For comparison only - Int_t label=track.GetLabel(); - if (label>0) ngood++; - - //cout << esdTrack->GetStatus() << " "; - esdTrack->UpdateTrackParams(&track,AliESDtrack::kITSrefit); - //cout << esdTrack->GetStatus() << endl; - ntrk++; - } - } - - Info("RefitInward","Refitted tracks: %d",ntrk); - if (ntrk) - Info("RefitInward","Good tracks/refitted: %f",Float_t(ngood)/ntrk); - - return 0; -} - -//________________________________________________________________________________ -Int_t AliITSUTrackerSA::LoadClusters(TTree *cluTree) { - //-------------------------------------------------------------------- - // This function reads the ITSU clusters from the tree, - // sort them, distribute over the internal tracker arrays, etc - //-------------------------------------------------------------------- - fITS->LoadClusters(cluTree); - fITS->ProcessClusters(); - // - for(int iL=0; iL<7; ++iL) { - fClustersTC[iL]=*fITS->GetLayerActive(iL)->GetClustersAddress(); - AliITSURecoLayer* lr = fITS->GetLayerActive(iL) ; // assign the layer which the cluster belongs to - for(int iC=0;iCGetEntriesFast();++iC) { - const AliITSUClusterPix *cl = (AliITSUClusterPix*)fClustersTC[iL]->At(iC); - float pos[3]; - cl->GetGlobalXYZ(pos); - float phi = TMath::PiOver2(); - if(Abs(pos[0])>1e-9) { - phi=TMath::ATan2(pos[1]-GetY(),pos[0]-GetX()); - if(phi<0.f) phi+=TMath::TwoPi(); - } else if(pos[1]<0.f) phi *= 3.f ; - AliITSURecoSens* sens = lr->GetSensorFromID(cl->GetVolumeId()); - const float alpha = sens->GetPhiTF(); - const float cov[3]={cl->GetSigmaZ2(),cl->GetSigmaYZ(),cl->GetSigmaY2()}; - - fLayer[iL].AddPoint(pos,cov,phi,alpha); - for ( int i=0 ; i<3; ++i ) { - fLayer[iL].Points.back().Label[i] = (cl->GetLabel(i)<0) ? -1 : cl->GetLabel(i); - } - } - - fLayer[iL].Sort(); - - } - return 0; -} - -//________________________________________________________________________________ -void AliITSUTrackerSA::UnloadClusters() { - //-------------------------------------------------------------------- - // This function unloads ITSU clusters from the RAM - //-------------------------------------------------------------------- - for(int i=0;i<7;++i) - fLayer[i].Clear(); - for(int i=0;i<5;++i) - fCells[i].clear(); - -} - -//________________________________________________________________________________ -void AliITSUTrackerSA::CellularAutomaton(AliESDEvent *event) { - - // Here it's implemented the Cellular Automaton routine - // Firstly the level of each doublet is set according to the level of - // the neighbour doublets. - // Doublet are considered to be neighbour if they share one point and the - // phi and theta direction difference of the two is below a cut value. - Int_t ntrk=0,ngood=0; - - for( int iL = 1; iL < 5; ++iL ) { - for( size_t iC1 = 0; iC1 < fCells[iL-1].size(); ++iC1 ) { - for( size_t iC2 = 0; iC2 < fCells[iL].size(); ++iC2 ) { - if( fCells[iL-1][iC1].Points[1]==fCells[iL][iC2].Points[0] && - fCells[iL-1][iC1].Points[2]==fCells[iL][iC2].Points[1] && - fCells[iL-1][iC1].Level >= fCells[iL][iC2].Level - 1 ) { - // The implementation of the curvature based matching has to be studied. - fCells[iL][iC2].Level = fCells[iL-1][iC1].Level+1; - fCells[iL][iC2].Neighbours.push_back(iC1); - } - } - } - } - - for (int level = 5; level > 0; --level ) { - vector roads; - roads.reserve(100); // should reserve() be based on number of clusters on outermost layer? - for (int iCL=4; iCL >= level-1; --iCL ) { - for (size_t iCell = 0; iCell < fCells[iCL].size(); ++iCell) { - if ( fCells[iCL][iCell].Level != level ) continue; - roads.push_back( Road(iCL,iCell) ); - for( size_t iN=0;iN 0 ) roads.push_back(roads.back()); - CandidatesTreeTraversal(roads,neigh,currD); - } - fCells[iCL][iCell].Level = -1; - } - } - - // for(size_t iR=0; iR-1) { - // roads[iR].Label = -1; - // if(fCells[i][roads[iR][i]].Label==-1) roads[iR].Label--; - // } - - // cout << fCells[i][roads[iR][i]].Label << " "; - // } - // } - // cout << " | " << roads[iR].Label << " | " << roads[iR].N << endl; - // } - vector candidates; - candidates.reserve(roads.size()); - - for (size_t iR=0; iR index; - index.reserve(candidates.size()); - for ( size_t i = 0; i < candidates.size(); ++i ) index.push_back(i); - Comparison comp(&candidates); - sort(index.begin(),index.end(),comp); - - for ( size_t cand = 0; cand < candidates.size(); ++cand ) { - const int ii = index[cand]; - - if ( candidates[ii].GetChi2() < 0. ) continue; - - // cout << candidates[ii].GetChi2() << " " << candidates[ii].GetNumberOfClusters() << " | " << candidates[ii].GetLabel() << " | "; - // for(int i=0;iGetLabel(0) << " "; - // } - // cout << endl; - - if( candidates[ii].GetChi2()/candidates[ii].GetNumberOfClusters() > fgkChi2Cut ) { - break; - } - bool goodTrack = true; - for ( int point = 0; point < candidates[ii].GetNumberOfClusters(); ++point ) { - int layer = (candidates[ii].GetClusterIndex(point)&0xf0000000)>>28; - int ind = (candidates[ii].GetClusterIndex(point)&0x0fffffff); - - if( (fLayer[ layer ].Points[ ind ].Used ) ) { - goodTrack = false; - } - - } - if(!goodTrack) { - continue; - } - for ( int point = 0; point < candidates[ii].GetNumberOfClusters(); ++point ) { - int layer = (candidates[ii].GetClusterIndex(point)&0xf0000000)>>28; - int ind = (candidates[ii].GetClusterIndex(point)&0x0fffffff); - fLayer[ layer ].Points[ ind ].Used = true; - } - - AliESDtrack outTrack; - CookLabel((AliKalmanTrack*)&candidates[ii],0.f); - ntrk++; - if(candidates[ii].GetChi2()>0) ngood++; - - // cout << candidates[ii].GetChi2() << " " << candidates[ii].GetNumberOfClusters() << " | " << candidates[ii].GetLabel() << " | "; - // for(int i=0;iGetLabel(0) << " "; - // } - // cout << endl; - - outTrack.UpdateTrackParams((AliKalmanTrack*)&candidates[ii],AliESDtrack::kITSin); - outTrack.SetLabel(candidates[ii].GetLabel()); - event->AddTrack(&outTrack); - } - } - Info("Clusters2Tracks","Reconstructed tracks: %d",ntrk); - if (ntrk) - Info("Clusters2Tracks","Good tracks/reconstructed: %f",Float_t(ngood)/ntrk); -} - - -//________________________________________________________________________________ -void AliITSUTrackerSA::CellsCreation(const int &cutLevel) { - // Make associations between two points on adjacent layers within an azimuthal window. - // Under consideration: - // - track parameter estimation using the primary vertex position - // To do: - // - last iteration - - float phiCut = 7.f; - if( cutLevel==0 ) phiCut = fPhiCut; - - // Doublets creation - vector doublets[6]; - for( int iL = 0 ; iL < 6 ; ++iL ) { - for ( int iC1 = 0 ; iC1 < fLayer[iL].N ; ++iC1 ) { - for ( int iC2 = 0; iC2 < fLayer[iL+1].N ; ++iC2 ) { - const float dPhi = Abs( fLayer[iL][iC1].Phi - fLayer[iL+1][iC2].Phi ); - if( dPhi < phiCut || Abs( dPhi - TMath::TwoPi() ) < phiCut) { - doublets[iL].push_back(Cell(iC1,iC2)); - if(Abs(fLayer[iL][iC1].XYZ[0]-fLayer[iL+1][iC2].XYZ[0])<1e-32) { - doublets[iL].back().Param[0] = 1e32; - } else { - doublets[iL].back().Param[0] = (fLayer[iL][iC1].XYZ[1]-fLayer[iL+1][iC2].XYZ[1])/(fLayer[iL][iC1].XYZ[0]-fLayer[iL+1][iC2].XYZ[0]); - } - const float r1 = Sqrt(fLayer[iL][iC1].XYZ[0] * fLayer[iL][iC1].XYZ[0] + fLayer[iL][iC1].XYZ[1] * fLayer[iL][iC1].XYZ[1]); - const float r2 = Sqrt(fLayer[iL+1][iC2].XYZ[0] * fLayer[iL+1][iC2].XYZ[0] + fLayer[iL+1][iC2].XYZ[1] * fLayer[iL+1][iC2].XYZ[1]); - doublets[iL].back().Param[1] = (fLayer[iL][iC1].XYZ[2]-fLayer[iL+1][iC2].XYZ[2])/(r1-r2); - doublets[iL].back().Label=-1; - for(int i=0;i<3;++i) { - for(int j=0;j<3;++j) { - if(fLayer[iL][iC1].Label[i]>-1&&fLayer[iL][iC1].Label[i]==fLayer[iL+1][iC2].Label[j]) - doublets[iL].back().Label = fLayer[iL][iC1].Label[i]; - } - } - } else if( fLayer[iL+1][iC2].Phi - fLayer[iL][iC1].Phi > phiCut ) break; - } - - } - } - - // Triplets creation - for( int iL = 5; iL > 0; --iL ) { - fCells[iL-1].clear(); - for ( size_t iD2 = 0; iD2 < doublets[iL].size(); ++iD2 ) { - for ( size_t iD1 = 0; iD1 < doublets[iL-1].size(); ++iD1 ) { - const int id1 = doublets[iL-1][iD1].Points[1]; - const int id2 = doublets[iL][iD2].Points[0]; - if ( id1 == id2 ) { - const int id3 = doublets[iL][iD2].Points[1]; - const float r3 = Sqrt( fLayer[iL+1][id3].XYZ[0] * fLayer[iL+1][id3].XYZ[0] + fLayer[iL+1][id3].XYZ[1] * fLayer[iL+1][id3].XYZ[1] ); - const float r2 = Sqrt( fLayer[iL][id2].XYZ[0] * fLayer[iL][id2].XYZ[0] + fLayer[iL][id2].XYZ[1] * fLayer[iL][id2].XYZ[1] ); - const float extrZ3 = doublets[iL-1][iD1].Param[1] * ( r3 - r2 ) + fLayer[iL][id2].XYZ[2] ; - const int iii = doublets[iL-1][iD1].Points[0]; - if ( Abs ( extrZ3 - fLayer[iL+1][id3].XYZ[2] ) < fZCut ) { - fCells[iL-1].push_back(Cell(doublets[iL-1][iD1].Points[0],id2,id3)); - fCells[iL-1].back().Param[0] = Curvature(fLayer[iL+1][id3].XYZ[0],fLayer[iL+1][id3].XYZ[1],fLayer[iL][id2].XYZ[0],fLayer[iL][id2].XYZ[1],fLayer[iL-1][iii].XYZ[0],fLayer[iL-1][iii].XYZ[1]); - fCells[iL-1].back().Param[1] = doublets[iL][iD2].Param[1]; - if(doublets[iL-1][iD1].Label==doublets[iL][iD2].Label&&doublets[iL][iD2].Label!=-1) - fCells[iL-1].back().Label=doublets[iL][iD2].Label; - else - fCells[iL-1].back().Label=-1; - } - } - } - } - } - -} - -//______________________________________________________________________________ -Bool_t AliITSUTrackerSA::InitTrackParams(AliITSUTrackCooked &track, int points[]) -{ - // Set the initial guess on track kinematics for propagation. - // Assume at least 3 points available - int lrOcc[AliITSUAux::kMaxLayers], nCl=0; - // - // we will need endpoints and middle layer - for (int i=fITS->GetNLayersActive()-1; i>=0; i--) { - if (points[i<<0x1]>-1) { - lrOcc[nCl++] = i; - track.SetClusterIndex(i,points[i<<0x1]); - } - } - - if (nCl<3) { - AliError(Form("Cannot estimate momentum of tracks with %d clusters",nCl)); - return kFALSE; - } - // - const int lr0 = lrOcc[0]; - const int lr1 = lrOcc[nCl/2]; - const int lr2 = lrOcc[nCl-1]; - // - const SpacePoint& cl0 = fLayer[lr0].Points[ points[lr0<<1] ]; - const SpacePoint& cl1 = fLayer[lr1].Points[ points[lr1<<1] ]; - const SpacePoint& cl2 = fLayer[lr2].Points[ points[lr2<<1] ]; - double cv = Curvature(cl0.XYZ[0],cl0.XYZ[1], cl1.XYZ[0],cl1.XYZ[1], cl2.XYZ[0],cl2.XYZ[1]); - - double tgl = (cl2.XYZ[2]-cl0.XYZ[2])/TMath::Sqrt((cl2.XYZ[0]-cl0.XYZ[0])*(cl2.XYZ[0]-cl0.XYZ[0])+(cl2.XYZ[1]-cl0.XYZ[1])*(cl2.XYZ[1]-cl0.XYZ[1])); - // - AliITSUClusterPix* clus = (AliITSUClusterPix*)fClustersTC[ lr0 ]->At( points[lr0<<1] ); - AliITSURecoLayer* lr = fITS->GetLayerActive(lr0); - AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId()); - double x = sens->GetXTF() + clus->GetX(); - double alp = sens->GetPhiTF(); - // printf("Alp: %f phi: %f\n",alp,phi); - double par[5] = {clus->GetY(),clus->GetZ(),0,tgl,cv}; - double cov[15] = { - 5*5, - 0, 5*5, - 0, 0 , 0.7*0.7, - 0, 0, 0, 0.7*0.7, - 0, 0, 0, 0, 10 - }; - track.Set(x,alp,par,cov); - - return kTRUE; -} - -//______________________________________________________________________________ -void AliITSUTrackerSA::CandidatesTreeTraversal(vector &candidates, const int &iD, const int &doubl) { - - if ( doubl < 0 ) return; - - candidates.back().AddElement(doubl,iD); - const int currentN = candidates.back().N; - for ( size_t iN = 0; iN < fCells[doubl][iD].Neighbours.size(); ++iN ) { - const int currD = doubl - 1 ; - const int neigh = fCells[doubl][iD].Neighbours[iN]; - - if ( iN > 0 ) { - candidates.push_back(static_cast(candidates.back())); - candidates.back().N = currentN; - } - - CandidatesTreeTraversal(candidates,neigh,currD); - } - - fCells[doubl][iD].Level = -1; - -} - -//______________________________________________________________________________ -Double_t AliITSUTrackerSA::RefitTrack(AliExternalTrackParam* trc, - Int_t clInfo[2*AliITSUAux::kMaxLayers], - Double_t rDest, Int_t stopCond) -{ - // refit track till radius rDest. - // if stopCond<0 : propagate till last cluster then stop - // if stopCond==0: propagate till last cluster then try to go till limiting rDest, don't mind if fail - // if stopCond>0 : rDest must be reached - // - // The clList should provide the indices of clusters at corresponding layer (as stored in the layer - // TClonesArray, with convention (allowing for up to 2 clusters per layer due to the overlaps): - // if there is a cluster on given layer I, then it should be stored at clInfo[2*I-1] - // if there is an additional cluster on this layer, it goes to clInfo[2*I] - // -1 means no cluster - // - double rCurr = Sqrt(trc->GetX()*trc->GetX() + trc->GetY()*trc->GetY()); - int dir,lrStart,lrStop; - // - dir = rCurrFindFirstLayerID(rCurr,dir); - lrStop = fITS->FindLastLayerID(rDest,dir); // lr id before which we have to stop - // - if (lrStop<0 || lrStart<0) AliFatal(Form("Failed to find start(%d) or last(%d) layers. " - "Track from %.3f to %.3f",lrStart,lrStop,rCurr,rDest)); - // - int nCl = 0; - for (int i=2*fITS->GetNLayersActive();i--;) {if (clInfo[i]<0) continue; nCl++;} - // - AliExternalTrackParam tmpTr(*trc); - double chi2 = 0; - int iclLr[2],nclLr; - int nclFit = 0; - // - int lrStop1 = lrStop+dir; - for (int ilr=lrStart;ilr!=lrStop1;ilr+=dir) { - AliITSURecoLayer* lr = fITS->GetLayer(ilr); - if ( dir*(rCurr-lr->GetR(dir))>0) continue; // this layer is already passed - int ilrA2,ilrA = lr->GetActiveID(); - // passive layer or active w/o hits will be traversed on the way to next cluster - if (!lr->IsActive() || clInfo[ilrA2=(ilrA<<1)]<0) continue; - // - // select the order in which possible 2 clusters (in case of the overlap) will be traversed and fitted - nclLr=0; - if (dir>0) { // clusters are stored in increasing radius order - iclLr[nclLr++]=clInfo[ilrA2++]; - if (clInfo[ilrA2]>=0) iclLr[nclLr++]=clInfo[ilrA2]; - } - else { - if ( clInfo[ilrA2+1]>=0 ) iclLr[nclLr++]=clInfo[ilrA2+1]; - iclLr[nclLr++]=clInfo[ilrA2]; - } - // - Bool_t transportedToLayer = kFALSE; - for (int icl=0;iclGetCluster(iclLr[icl]); - AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId()); - if (!tmpTr.Rotate(sens->GetPhiTF())) return -1; - // - double xClus = sens->GetXTF()+clus->GetX(); - if (!transportedToLayer) { - if (ilr!=lrStart && !TransportToLayerX(&tmpTr,lrStart,ilr,xClus)) return -1; // go to the entrance to the layer - lrStart = ilr; - transportedToLayer = kTRUE; - } - // - if (!PropagateSeed(&tmpTr,xClus,fCurrMass)) return -1; - // - Double_t p[2]={clus->GetY(), clus->GetZ()}; - Double_t cov[3]={clus->GetSigmaY2(), clus->GetSigmaYZ(), clus->GetSigmaZ2()}; - double chi2cl = tmpTr.GetPredictedChi2(p,cov); - chi2 += chi2cl; - // - if ( !tmpTr.Update(p,cov) ) return -1; - if (++nclFit==nCl && stopCond<0) { - *trc = tmpTr; - return chi2; // it was requested to not propagate after last update - } - } - // - } - // All clusters were succesfully fitted. Even if the track does not reach rDest, this is enough to validate it. - // Still, try to go as close as possible to rDest. - // - if (lrStart!=lrStop) { - if (!TransportToLayer(&tmpTr,lrStart,lrStop)) return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory - if (!GoToExitFromLayer(&tmpTr,fITS->GetLayer(lrStop),dir)) return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory - } - // go to the destination radius. Note that here we don't select direction to avoid precision problems - if (!tmpTr.GetXatLabR(rDest,rDest,GetBz(),0) || !PropagateSeed(&tmpTr,rDest,fCurrMass, 100, kFALSE)) { - return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory - } - *trc = tmpTr; - - return chi2; -} - -//______________________________________________________________________________ -Bool_t AliITSUTrackerSA::PropagateSeed(AliExternalTrackParam *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr) -{ - // propagate seed to given x applying material correction if requested - const Double_t kEpsilon = 1e-5; - Double_t xpos = seed->GetX(); - Int_t dir = (xpos0 && seed->IsStartedTimeIntegral(); - if (matCorr || updTime) seed->GetXYZ(xyz1); //starting global position - while ( (xToGo-xpos)*dir > kEpsilon){ - Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep); - Double_t x = xpos+step; - Double_t bz=GetBz(); // getting the local Bz - if (!seed->PropagateTo(x,bz)) return kFALSE; - double ds = 0; - if (matCorr || updTime) { - xyz0[0]=xyz1[0]; // global pos at the beginning of step - xyz0[1]=xyz1[1]; - xyz0[2]=xyz1[2]; - seed->GetXYZ(xyz1); // // global pos at the end of step - // - if (matCorr) { - Double_t xrho,xx0; - ds = GetMaterialBudget(xyz0,xyz1,xx0,xrho); - if (dir>0) xrho = -xrho; // outward should be negative - if (!seed->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE; - } - else { // matCorr is not requested but time integral is - double d0 = xyz1[0]-xyz0[0]; - double d1 = xyz1[1]-xyz0[1]; - double d2 = xyz1[2]-xyz0[2]; - ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2); - } - } - if (updTime) seed->AddTimeStep(ds); - // - xpos = seed->GetX(); - } - return kTRUE; -} - -//_________________________________________________________________________ -Bool_t AliITSUTrackerSA::TransportToLayer(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t rLim) -{ - // transport track from layerFrom to the entrance of layerTo or to rLim (if>0), wathever is closer - // - if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo)); - // - int dir = lTo > lFrom ? 1:-1; - AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested - Bool_t checkFirst = kTRUE; - Bool_t limReached = kFALSE; - while(lFrom!=lTo) { - if (lrFr) { - if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer - checkFirst = kFALSE; - } - AliITSURecoLayer* lrTo = fITS->GetLayer( (lFrom+=dir) ); - if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom)); - // - // go the entrance of the layer, assuming no materials in between - double xToGo = lrTo->GetR(-dir); - if (rLim>0) { - if (dir>0) { - if (rLimxToGo) {xToGo = rLim; limReached = kTRUE;} - } - } - // double xts = xToGo; - if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) { - // printf("FailHere1: %f %f %d\n",xts,xToGo,dir); - // seed->Print("etp"); - return kFALSE; - } - if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) { - //printf("FailHere2: %f %f %d\n",xts,xToGo,dir); - //seed->Print("etp"); - return kFALSE; - } - lrFr = lrTo; - if (limReached) break; - } - return kTRUE; - // -} - -//_________________________________________________________________________ -Bool_t AliITSUTrackerSA::TransportToLayerX(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t xStop) -{ - // transport track from layerFrom to the entrance of layerTo but do not pass control parameter X - // - if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo)); - // - int dir = lTo > lFrom ? 1:-1; - AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested - Bool_t checkFirst = kTRUE; - while(lFrom!=lTo) { - if (lrFr) { - if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer - checkFirst = kFALSE; - } - AliITSURecoLayer* lrTo = fITS->GetLayer( (lFrom+=dir) ); - if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom)); - // - // go the entrance of the layer, assuming no materials in between - double xToGo = lrTo->GetR(-dir); // R of the entrance to layer - // - // double xts = xToGo; - if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) { - // printf("FailHere1: %f %f %d\n",xts,xToGo,dir); - // seed->Print("etp"); - return kFALSE; - } - if ( (dir>0&&xToGo>xStop) || (dir<0&&xToGoGetR(-dir), xToGo)); -#endif - if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) { - //printf("FailHere2: %f %f %d\n",xts,xToGo,dir); - //seed->Print("etp"); - return kFALSE; - } - lrFr = lrTo; - } - return kTRUE; - // -} - -//_________________________________________________________________________ -Bool_t AliITSUTrackerSA::GoToExitFromLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check) -{ - // go to the exit from lr in direction dir, applying material corrections in steps specific for this layer - // If check is requested, do this only provided the track has not exited the layer already - double xToGo = lr->GetR(dir); - if (check) { // do we need to track till the surface of the current layer ? - double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius - if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // on the surface or outside of the layer - else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // on the surface or outside of the layer - } - if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE; - // go via layer to its boundary, applying material correction. - if (!PropagateSeed(seed,xToGo,fCurrMass, lr->GetMaxStep())) return kFALSE; - // - return kTRUE; - // -} - -//_________________________________________________________________________ -Bool_t AliITSUTrackerSA::GoToEntranceToLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check) -{ - // go to the entrance of lr in direction dir, w/o applying material corrections. - // If check is requested, do this only provided the track did not reach the layer already - double xToGo = lr->GetR(-dir); - if (check) { // do we need to track till the surface of the current layer ? - double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius - if (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // already passed - else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // already passed - } - if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE; - // go via layer to its boundary, applying material correction. - if (!PropagateSeed(seed,xToGo,fCurrMass, 100, kFALSE)) return kFALSE; - return kTRUE; - // -} - -//____________________________________________________ -Double_t AliITSUTrackerSA::GetMaterialBudget(const double* pnt0,const double* pnt1, double& x2x0, double& rhol) const -{ - double par[7]; - if (fUseMatLUT && fMatLUT) { - double d = fMatLUT->GetMatBudget(pnt0,pnt1,par); - x2x0 = par[AliITSUMatLUT::kParX2X0]; - rhol = par[AliITSUMatLUT::kParRhoL]; - return d; - } - else { - MeanMaterialBudget(pnt0,pnt1,par); - x2x0 = par[1]; - rhol = par[0]*par[4]; - return par[4]; - } -} - -//____________________________________________________________________ -Double_t AliITSUTrackerSA::Curvature(Double_t x1,Double_t y1,Double_t x2,Double_t y2,Double_t x3,Double_t y3) { - - //calculates the curvature of track - Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1); - if(den*den<1e-32) return 0.; - Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den; - if((y2-y1)*(y2-y1)<1e-32) return 0.; - Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1); - Double_t c = -x1*x1-y1*y1-a*x1-b*y1; - Double_t xc= -a/2.; - - if((a*a+b*b-4*c)<0) return 0.; - Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.; - if(rad*rad<1e-32) return 1e16; - - if((x1>0 && y1>0 && x10 && x10 && y1<0 && x1 -#include -#include "AliITSUTrackerSAaux.h" // Structs and other stuff -#include "AliITSUMatLUT.h" -#include "AliITSUAux.h" -#include "AliExternalTrackParam.h" - -class AliITSUReconstructor; -class AliITSURecoDet; -class AliITSURecoLayer; - -class TTree; -class AliCluster; -class AliESDEvent; - -using std::vector; - -//------------------------------------------------------------------------- -class AliITSUTrackerSA : public AliTracker { -public: - AliITSUTrackerSA(AliITSUReconstructor* rec=0); - virtual ~AliITSUTrackerSA(); - - // These functions must be implemented - Int_t Clusters2Tracks(AliESDEvent *event); - Int_t PropagateBack(AliESDEvent *event); - Int_t RefitInward(AliESDEvent *event); - Int_t LoadClusters(TTree *ct); - void UnloadClusters(); - - inline AliCluster *GetCluster(Int_t index) const { - const Int_t l=(index & 0xf0000000) >> 28; - const Int_t c=(index & 0x0fffffff); - return (AliCluster*)fClustersTC[l]->At(c) ; - } - - // Possibly, other public functions - void Init(AliITSUReconstructor* rec); - Double_t RefitTrack(AliExternalTrackParam* trc, Int_t clInfo[2*AliITSUAux::kMaxLayers], Double_t rDest, Int_t stopCond); - Bool_t PropagateSeed(AliExternalTrackParam *seed, Double_t xToGo, Double_t mass, Double_t maxStep=1.0, Bool_t matCorr=kTRUE); - Double_t GetMaterialBudget(const double* pnt0, const double* pnt1, double& x2x0, double& rhol) const; - Bool_t GoToEntranceToLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check=kFALSE); - Bool_t GoToExitFromLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check=kTRUE); - Bool_t TransportToLayerX(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t xStop); - Bool_t TransportToLayer(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t rLim=-1); - - void SetChi2Cut(float cut) { fChi2Cut=cut; } - void SetRPhiCut(float cut) { fRPhiCut=cut; } - void SetPhiCut(float cut) { fPhiCut=cut; } - void SetZCut(float cut) { fZCut=cut; } - - - // -protected: - AliITSUTrackerSA(const AliITSUTrackerSA&); - - void CellsCreation(const int &cutLevel); - void CellularAutomaton(AliESDEvent *ev); - // void MakeTriplets(); - void CandidatesTreeTraversal( vector &vec, const int &iD, const int &doubl); - Bool_t InitTrackParams(AliITSUTrackCooked &track, int points[]); - void GlobalFit(); - void ChiSquareSelection(); - void MergeTracks( vector &vec, bool flags[] ); - // Other protected functions - // (Sorting, labeling, calculations of "roads", etc) - static Double_t Curvature(Double_t x1,Double_t y1,Double_t x2,Double_t y2,Double_t x3,Double_t y3); - - - -private: - AliITSUTrackerSA &operator=(const AliITSUTrackerSA &tr); - void SetLabel(AliITSUTrackCooked &t, Float_t wrong); - - // Data members - - // classes for interfacing the geometry, materials etc. - AliITSUReconstructor* fReconstructor; // ITS global reconstructor - AliITSURecoDet* fITS; // interface to ITS, borrowed from reconstructor - AliITSUMatLUT* fMatLUT; // material lookup table - Bool_t fUseMatLUT; //! use material lookup table rather than TGeo - Double_t fCurrMass; // assumption about particle mass - - - // Internal tracker arrays, layers, modules, etc - Layer fLayer[7]; - TClonesArray *fClustersTC[7]; - vector fCells[5]; - Float_t fChi2Cut; - Float_t fPhiCut; - Float_t fRPhiCut; - Float_t fZCut; - - // - static const Double_t fgkToler; // tracking tolerance - static const Double_t fgkChi2Cut; // chi2 cut during track merging - // - ClassDef(AliITSUTrackerSA,1) //ITSU stand-alone tracker -}; - -#endif diff --git a/ITS/UPGRADE/AliITSUTrackerSAaux.h b/ITS/UPGRADE/AliITSUTrackerSAaux.h deleted file mode 100644 index 11e87421830..00000000000 --- a/ITS/UPGRADE/AliITSUTrackerSAaux.h +++ /dev/null @@ -1,268 +0,0 @@ -#ifndef ALIITSUTRACKERSAAUX_H -#define ALIITSUTRACKERSAAUX_H - -#include -using std::vector; -#include -using std::sort; -#include "AliExternalTrackParam.h" -#include "AliITSUTrackCooked.h" -#include "AliITSUAux.h" - -#include - -#include -using std::cout; -using std::endl; - -//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: -template -struct Comparison { //Adapted from TMath ROOT code - Comparison(vector *d) : fData(d) {} - - bool operator()(int i1, int i2) { - return fData->at(i1).CmpVar() < fData->at(i2).CmpVar(); - } - vector *fData; -}; - -//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: -struct Cell { - - Cell() : Level(1), - Param(), - Points(), - N(0), - Neighbours(), - Label(-1){ - // Inlined standard constructor - } ; - - Cell(int i1, int i2) : Level(1), - Param(), - Points(), - N(2), - Neighbours(), - Label(-1){ - // Inlined standard constructor - Points[0] = i1; - Points[1] = i2; - } ; - - Cell(int i1, int i2, int i3) : Level(1), - Param(), - Points(), - N(3), - Neighbours(), - Label(-1) { - - // Inlined standard constructor - Points[0] = i1; - Points[1] = i2; - Points[2] = i3; - } ; - - int Level; - float Param[3]; - int Points[3]; - int N; - vector Neighbours; - int Label; -}; - -//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: -struct SpacePoint { - SpacePoint(const float *xyz, const float &alpha) : XYZ(), - Cov(), - Phi(), - Alpha(alpha), - Label(), - Used(false) - { - XYZ[0] = xyz[0]; - XYZ[1] = xyz[1]; - XYZ[2] = xyz[2]; - Cov[0]=Cov[1]=Cov[2]=99999.f; - }; - - SpacePoint(const SpacePoint &sp) : XYZ(), - Cov(), - Phi(sp.Phi), - Alpha(sp.Alpha), - Label(), - Used(sp.Used) - { - for(int i=0; i<3; ++i) { - XYZ[i]=sp.XYZ[i]; - Label[i]=sp.Label[i]; - Cov[i]=sp.Cov[i]; - } - }; - - float CmpVar() const { return Phi; } - - float XYZ[3]; - float Cov[3]; - float Phi; - float Alpha; - int Label[3]; - bool Used; -}; - -//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: -struct Layer { - Layer() : Index(), - N(0), - Points() {}; - - int operator()(const int &i) { return Index[i]; } // Just for compatibility with old code - SpacePoint& operator[](const int &i) { return Points[Index[i]]; } - ~Layer() { Clear(); } - - void AddPoint(const float xyz[3], const float cov[3], const float &phi, const float &alpha) { - Index.push_back(N++); - Points.push_back(SpacePoint(xyz,alpha)); - Points.back().Cov[0] = cov[0]; - Points.back().Cov[1] = cov[1]; - Points.back().Cov[2] = cov[2]; - Points.back().Phi = phi; - } - - void Sort() { - Comparison cmp(&Points); - sort(Index.begin(),Index.end(),cmp); - } - - void Clear() { - Index.clear(); - Points.clear(); - N=0; - } - - vector Index; - int N; - vector Points; -}; - -//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: -struct Road { - - Road() : Elements(), N(0), Label(-1) { - ResetElements(); - } - - Road(int layer, int id) : Elements(), N(1), Label(-1) { - ResetElements(); - Elements[layer] = id; - } - - Road(const Road& copy) : Elements(), N(copy.N), Label(copy.Label) { - for ( int i=0; i<5; ++i ) Elements[i] = copy.Elements[i]; - } - - int &operator[] (const int &i) { - return Elements[i]; - } - - void ResetElements() { - for ( int i=0; i<5; ++i ) - Elements[i] = -1; - } - - void AddElement(int i, int el) { - ++N; - Elements[i] = el; - } - - int Elements[5]; - int N; - int Label; -}; - -// class AliITSUTrackCA : public AliITSUTrackCooked { -// public: -// AliITSUTrackCA() : AliITSUTrackCooked() {} -// AliITSUTrackCA(const AliITSUTrackCA &t) : AliITSUTrackCooked((AliITSUTrackCooked)t) {} -// AliITSUTrackCA(const AliESDtrack &t) : AliITSUTrackCooked(t) {} -// AliITSUTrackCA &operator=(const AliITSUTrackCooked &t); -// virtual ~AliITSUTrackCA(); -// void SetChi2(Double_t chi2) { fChi2=chi2; } -// ClassDef(AliITSUTrackCA,1) -// }; - -//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: -template<> -struct Comparison { - Comparison(vector *d) : fData(d) {} - - bool operator()(int i1, int i2) { - return (fData->at(i1).GetChi2() < fData->at(i2).GetChi2()); - } - vector *fData; -}; - -// //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: -// class Candidate : public AliKalmanTrack { -// public : - -// Candidate() : AliKalmanTrack(), -// fChi2( 0. ), -// fPoints(), -// fNPoints(0), -// fInnermostLayer(-1), -// fOutermostLayer(-1) -// { -// for ( unsigned int i = 0; i < 2*AliITSUAux::kMaxLayers; ++i ) fPoints[i]=-1; -// }; - -// Candidate(const Candidate ©) : AliExternalTrackParam(), -// fChi2(copy.fChi2), -// fPoints(), -// fNPoints(copy.fNPoints), -// fInnermostLayer(copy.fInnermostLayer), -// fOutermostLayer(copy.fOutermostLayer) -// { -// for ( unsigned int i = 0; i < 2*AliITSUAux::kMaxLayers; ++i ) fPoints[i]=copy.fPoints[i]; -// }; - -// Candidate(int points[7]) : AliExternalTrackParam(), -// fChi2( 0. ), -// fPoints(), -// fNPoints( 0 ), -// fInnermostLayer( -1 ), -// fOutermostLayer( -1 ) -// { -// bool outer=false; -// ResetPoints(); -// for ( int i=6; i>=0; --i ) { -// if (points[i]!=-1) { -// if (! outer ) { -// outer = true; -// fOutermostLayer = points[i]; -// } -// fInnermostLayer = points[i]; -// ++fNPoints; -// } -// fPoints[i<<0x1] = points[i]; -// } -// } - -// int GetClusterIndex(const int &i) const { -// return ((fInnermostLayer+i)<<28)+fPoints[(fInnermostLayer+i)<<0x1]; -// } - -// Double_t CmpVar() const { return fChi2; } - -// void ResetPoints() { for(unsigned int i = 0; i < 2*AliITSUAux::kMaxLayers; ++i) fPoints[i]=-1; } -// // -// // Double_t fChi2; -// // Double_t fFakeRatio; -// Int_t fPoints[2*AliITSUAux::kMaxLayers]; -// // Int_t fNPoints; -// Int_t fLabel; -// Int_t fInnermostLayer; -// Int_t fOutermostLayer; - -// }; - -#endif diff --git a/ITS/UPGRADE/CMakelibITSUpgradeRec.pkg b/ITS/UPGRADE/CMakelibITSUpgradeRec.pkg index 03bd904add4..a84ba8585fe 100644 --- a/ITS/UPGRADE/CMakelibITSUpgradeRec.pkg +++ b/ITS/UPGRADE/CMakelibITSUpgradeRec.pkg @@ -37,7 +37,8 @@ set ( SRCS AliITSUClusterPix.cxx AliITSUSeed.cxx AliITSUTrackerGlo.cxx - AliITSUTrackerSA.cxx + AliITSUCATracker.cxx + AliITSUCATrackingStation.cxx AliITSUTrackerCooked.cxx AliITSUTrackCooked.cxx AliITSUTrackCond.cxx @@ -61,8 +62,6 @@ string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" ) set ( DHDR ITSUpgradeRecLinkDef.h) -#set_source_files_properties( AliITSUTrackerSA.cxx PROPERTIES COMPILE_FLAGS "-O3 -ftree-vectorize -ffast-math -std=gnu++11" ) - set ( EINCLUDE TPC RAW ITS ITS/UPGRADE ITS/UPGRADE/v0 STEER/STEER STEER/ESD STEER/STEERBase) # set ( EXPORT AliITStrackV2.h AliITSVertexer.h AliITSRecoParam.h) diff --git a/ITS/UPGRADE/ITSUpgradeRecLinkDef.h b/ITS/UPGRADE/ITSUpgradeRecLinkDef.h index 4ce519cb261..640fbf81d0d 100644 --- a/ITS/UPGRADE/ITSUpgradeRecLinkDef.h +++ b/ITS/UPGRADE/ITSUpgradeRecLinkDef.h @@ -23,7 +23,7 @@ #pragma link C++ class AliITSUClusterPix+; #pragma link C++ class AliITSUSeed+; #pragma link C++ class AliITSUTrackerGlo+; -#pragma link C++ class AliITSUTrackerSA+; +#pragma link C++ class AliITSUCATracker+; #pragma link C++ class AliITSUTrackerCooked+; #pragma link C++ class AliITSUTrackCooked+; #pragma link C++ class AliITSUTrackCond+; diff --git a/ITS/UPGRADE/macros/testTrackerCA.C b/ITS/UPGRADE/macros/testTrackerCA.C index 51277191709..5f0869cf269 100644 --- a/ITS/UPGRADE/macros/testTrackerCA.C +++ b/ITS/UPGRADE/macros/testTrackerCA.C @@ -6,6 +6,8 @@ #include #include #include +#include +#include #include "AliTracker.h" #include "AliGeomManager.h" @@ -19,7 +21,7 @@ #include "AliESDEvent.h" #include "AliITSURecoParam.h" #include "AliITSUReconstructor.h" -#include "AliITSUTrackerSA.h" +#include "AliITSUCATracker.h" #endif extern TSystem *gSystem; @@ -63,7 +65,10 @@ void testTrackerCA() { rec->SetRecoParam(par); // rec->Init(); - AliITSUTrackerSA *tracker = new AliITSUTrackerSA(); + AliITSUCATracker *tracker = new AliITSUCATracker(); + tracker->SetPhiCut(0.15f); + //tracker->SetThetaBin(0.03f); + tracker->SetChi2Cut(3000); tracker->Init(rec); @@ -77,11 +82,14 @@ void testTrackerCA() { AliRunLoader *rl = AliRunLoader::Open("galice.root","something"); rl->LoadHeader(); - + + TH1F fGHisto("ghisto","Good tracks #chi^{2};#chi^{2};Entries",100,0,400); + TH1F fFHisto("fhisto","Fake tracks #chi^{2};#chi^{2};Entries",100,0,400); TStopwatch timer; - Int_t nEvents=1;//rl->GetNumberOfEvents(); + Int_t nEvents=rl->GetNumberOfEvents(); for (Int_t i=0; iGetEvent(i); @@ -91,8 +99,10 @@ void testTrackerCA() { TTree *cTree=(TTree *)clsFile->Get(Form("Event%d/TreeR",i)); tracker->LoadClusters(cTree); tracker->Clusters2Tracks(esd); - //tracker->PropagateBack(esd); + tracker->PropagateBack(esd); tracker->RefitInward(esd); + // fGHisto.Add(tracker->GetGHisto()); + // fFHisto.Add(tracker->GetFHisto()); tracker->UnloadClusters(); Int_t n=esd->GetNumberOfTracks(); @@ -109,13 +119,88 @@ void testTrackerCA() { delete vtx; } timer.Stop(); timer.Print(); - + #ifdef _TUNING_ + TCanvas *cv1 = new TCanvas("chi2","chi2",1400,1000); + cv1->Divide(2,2); + for(int i=0;i<4;i++) { + cv1->cd(i+1); + cv1->cd(i+1)->SetLogy(); + tracker->fGoodCombChi2[i]->Draw(); + tracker->fFakeCombChi2[i]->SetLineColor(kRed); + tracker->fFakeCombChi2[i]->Draw("same"); + } + TCanvas *ccv1 = new TCanvas("deltan","DeltaN",1400,1000); + ccv1->Divide(2,2); + for(int i=0;i<4;i++) { + ccv1->cd(i+1); + ccv1->cd(i+1)->SetLogy(); + tracker->fGoodCombN[i]->Draw(); + tracker->fFakeCombN[i]->SetLineColor(kRed); + tracker->fFakeCombN[i]->Draw("same"); + } + TCanvas *cv2 = new TCanvas("chi2tracks","chi2tracks"); + cv2->cd(); + cv2->cd()->SetLogy(); + tracker->fGoodCombChi2[4]->Draw(); + tracker->fFakeCombChi2[4]->SetLineColor(kRed); + tracker->fFakeCombChi2[4]->Draw("same"); + TCanvas *cv3 = new TCanvas(); + cv3->cd(); + tracker->fTanF->SetLineColor(kRed); + tracker->fTanF->Draw(""); + tracker->fTan->Draw("same"); + TCanvas *cv33 = new TCanvas(); + cv33->cd(); + tracker->fPhiF->SetLineColor(kRed); + tracker->fPhiF->Draw(""); + tracker->fPhi->Draw("same"); + // TCanvas *cv4 = new TCanvas(); + // cv4->cd(); + // tracker->fNEntries->Draw(); + for (int i = 0; i < 6; ++i) + { + TCanvas *cv = new TCanvas(Form("Doublets%i",i),Form("Doublets%i",i),1400,500); + cv->Divide(2); + cv->cd(1); + tracker->fGDZ[i]->Draw(); + tracker->fFDZ[i]->SetLineColor(kRed); + tracker->fFDZ[i]->Draw("same"); + cv->cd(2); + tracker->fGDXY[i]->Draw(); + tracker->fFDXY[i]->SetLineColor(kRed); + tracker->fFDXY[i]->Draw("same"); + } + for (int i = 0; i < 5; ++i) + { + TCanvas *cv = new TCanvas(Form("Cells%i",i),Form("Cells%i",i),1400,500); + cv->Divide(2); + cv->cd(1); + tracker->fGDCAZ[i]->Draw(); + tracker->fFDCAZ[i]->SetLineColor(kRed); + tracker->fFDCAZ[i]->Draw("same"); + cv->cd(2); + tracker->fGDCAXY[i]->Draw(); + tracker->fFDCAXY[i]->SetLineColor(kRed); + tracker->fFDCAXY[i]->Draw("same"); + } + // TCanvas *ccv4 = new TCanvas("iparams","Impact Parameters",1000,700); + // ccv4->Divide(2); + // ccv4->cd(1); + // tracker->fIPGoodXY->Draw(); + // tracker->fIPFakeXY->SetLineColor(kRed); + // tracker->fIPFakeXY->Draw("same"); + // ccv4->cd(2); + // tracker->fIPGoodZ->Draw(); + // tracker->fIPFakeZ->SetLineColor(kRed); + // tracker->fIPFakeZ->Draw("same"); + #endif delete tracker; //delete clsFile; esdFile->cd(); esdTree->Write(); delete esd; delete esdFile; + } const AliESDVertex *SetMCvertex(const AliRunLoader *rl, AliTracker *tracker) { -- 2.39.3