Big improvement in terms of CPU speed.
namespace AliITSUAux {
- void BringTo02Pi(double &phi);
+ template<typename F>
+ 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);
}
//_________________________________________________________________________________
-inline void AliITSUAux::BringTo02Pi(double &phi) {
+template<typename F>
+inline void AliITSUAux::BringTo02Pi(F &phi) {
// bring phi to 0-2pi range
if (phi<0) phi+=TwoPi(); else if (phi>TwoPi()) phi-=TwoPi();
}
--- /dev/null
+#ifndef ALIITSUCACELL_H
+#define ALIITSUCACELL_H
+
+#include <vector>
+#include <iostream>
+
+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<int> 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
--- /dev/null
+//-*- 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 <maximiliano.puccio@cern.ch> *
+// 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 <algorithm>
+#include <TBranch.h>
+#include <TMath.h>
+#include <TTree.h>
+#include <Riostream.h>
+#include "AliLog.h"
+#include "AliESDEvent.h"
+#include "AliITSUClusterPix.h"
+#include "AliITSUCATracker.h"
+#include "AliITSUReconstructor.h"
+#include "AliITSURecoDet.h"
+#include "AliESDtrack.h"
+#include <cassert>
+
+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 && x1<xc)) rad*=-1;
+ // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
+
+ return 1.f/rad;
+
+}
+
+//__________________________________________________________________________________________________
+static inline float TanLambda(float x1, float y1, float x2, float y2, float z1, float z2)
+{
+ //
+ // Initial approximation of the tangent of the track dip angle
+ //
+ return (z1 - z2) * invsqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
+}
+
+//__________________________________________________________________________________________________
+//static inline float XCenterOfCurvature(float x1, float y1, float x2, float y2, float x3, float y3)
+//{
+// //
+// // Initial approximation of the x-coordinate of the center of curvature
+// //
+//
+// const float k1 = (y2 - y1) / (x2 - x1), k2 = (y3 - y2) / (x3 - x2);
+// return TMath::Abs(k2 - k1) > 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<AliITSUCARoad> &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<int> 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<int> 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<AliITSUCARoad> 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; i<n; i++) {
+ AliESDtrack *esdTrack=event->GetTrack(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 = rCurr<rDest ? 1 : -1;
+ lrStart = fITS->FindFirstLayerID(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;icl<nclLr;icl++) {
+ AliITSUClusterPix* clus = fLayer[ilrA].GetClusterUnSorted(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 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;
+}
+
--- /dev/null
+//-*- 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 <TClonesArray.h>
+#include <vector>
+#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<AliITSUCARoad> &roads, const int &iD, const int &doubl);
+ Bool_t InitTrackParams(AliITSUTrackCooked &track, int points[]);
+ void FindTracksCA(int iteration);
+ void GlobalFit();
+ void MergeTracks( vector<AliITSUTrackCooked> &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<bool> fUsedClusters[7];
+ Float_t fChi2Cut;
+ Float_t fPhiCut;
+ Float_t fZCut;
+ vector<Doublets> fDoublets[6];
+ vector<AliITSUCACell> 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
--- /dev/null
+//-*- 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 <maximiliano.puccio@cern.ch> *
+// 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 <TMath.h>
+#include "AliITSUGeomTGeo.h"
+#include "AliVertex.h"
+#include <TRandom.h>
+#include <TStopwatch.h>
+#include <TString.h>
+#include "AliITSUAux.h"
+#include "AliITSURecoSens.h"
+#include <Riostream.h>
+
+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<fNOccBins;i++) {
+ printf("%4d %5d(phi/z: %03d/%03d) -> %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;
+}
--- /dev/null
+//-*- 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 <algorithm>
+#include <vector>
+#include <TObject.h>
+#include <TClonesArray.h>
+#include "AliITSURecoLayer.h"
+#include "AliITSUClusterPix.h"
+#include <assert.h>
+
+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<int> fIndex;
+ std::vector<int> fFoundBins; // occupied bins satisfying to query
+ std::vector<ClsInfo_t> fSortedClInfo; // processed cluster info
+ std::vector<ITSDetInfo_t> 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
+++ /dev/null
-//--------------------------------------------------------------------------------
-// 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 <TBranch.h>
-#include <TMath.h>
-using TMath::Abs;
-using TMath::Sort;
-using TMath::Sqrt;
-#include <TTree.h>
-#include <algorithm>
-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 <Riostream.h>
-
-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; 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 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; 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 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;iC<fClustersTC[iL]->GetEntriesFast();++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<Road> 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<fCells[iCL][iCell].Neighbours.size(); ++iN ) {
- const int currD = iCL - 1;
- const int neigh = fCells[iCL][iCell].Neighbours[iN];
- if( iN > 0 ) roads.push_back(roads.back());
- CandidatesTreeTraversal(roads,neigh,currD);
- }
- fCells[iCL][iCell].Level = -1;
- }
- }
-
- // for(size_t iR=0; iR<roads.size(); ++iR) {
- // cout << "ROAD " << iR << " | ";
- // for(int i=0;i<5;++i) {
- // if(roads[iR][i]<0) continue;
- // else {
- // if(roads[iR].Label==-1){
- // roads[iR].Label = fCells[i][roads[iR][i]].Label;
- // if(roads[iR].Label==-1) roads[iR].Label--;
- // }
- // if (fCells[i][roads[iR][i]].Label!=roads[iR].Label&&roads[iR].Label>-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<AliITSUTrackCooked> candidates;
- candidates.reserve(roads.size());
-
- for (size_t iR=0; iR<roads.size(); ++iR) {
- if(roads[iR].N != level) {
- continue;
- }
-
- int points[2*AliITSUAux::kMaxLayers];
- for(unsigned int i=0;i<2*AliITSUAux::kMaxLayers;++i) points[i] = -1;
- for(int i=0;i<5;++i) {
- if(roads[iR].Elements[i]<0) continue;
- points[( i )<<1]=fLayer[ i ](fCells[i][roads[iR].Elements[i]].Points[0]);
- points[(i+1)<<1]=fLayer[i+1](fCells[i][roads[iR].Elements[i]].Points[1]);
- points[(i+2)<<1]=fLayer[i+2](fCells[i][roads[iR].Elements[i]].Points[2]);
- }
-
- candidates.push_back(AliITSUTrackCooked());
-
- InitTrackParams(candidates.back(),points);
- const double chi2 = RefitTrack( (AliExternalTrackParam*)&candidates.back(), points, 0. ,-1);
-
- if ( chi2 < 0. ) {
- // cout << "FAIL: " << chi2 << endl;
- // for(unsigned int i=0;i<2*AliITSUAux::kMaxLayers;++i)
- // cout << points[i] << " ";
- // cout << endl;
- candidates.back().SetChi2( 1e27 );
- } else candidates.back().SetChi2( chi2 );
- candidates.back().SetLabel(roads[iR].Label);
- }
-
- vector<int> index;
- index.reserve(candidates.size());
- for ( size_t i = 0; i < candidates.size(); ++i ) index.push_back(i);
- Comparison<AliITSUTrackCooked> 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;i<candidates[ii].GetNumberOfClusters();++i) {
- // cout<< GetCluster(candidates[ii].GetClusterIndex(i))->GetLabel(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;i<candidates[ii].GetNumberOfClusters();++i) {
- // cout<< GetCluster(candidates[ii].GetClusterIndex(i))->GetLabel(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<Cell> 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<Road> &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<Road>(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 = rCurr<rDest ? 1 : -1;
- lrStart = fITS->FindFirstLayerID(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;icl<nclLr;icl++) {
- AliITSUClusterPix* clus = (AliITSUClusterPix*)lr->GetCluster(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 = (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 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 (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 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&&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 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 && x1<xc)) rad*=-1;
- if((x1<0 && y1>0 && x1<xc)) rad*=-1;
- // if((x1<0 && y1<0 && x1<xc)) rad*=-1;
- // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
-
- return 1/rad;
-
-}
-
+++ /dev/null
-#ifndef ALIITSUTRACKERSA_H
-#define ALIITSUTRACKERSA_H
-
-//-------------------------------------------------------------------------
-// The stand-alone ITSU tracker
-// It reads AliITSUClusterPix clusters and writes the tracks to the ESD
-//-------------------------------------------------------------------------
-
-#include "AliTracker.h"
-#include "AliITSUGeomTGeo.h"
-#include <TClonesArray.h>
-#include <vector>
-#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<Road> &vec, const int &iD, const int &doubl);
- Bool_t InitTrackParams(AliITSUTrackCooked &track, int points[]);
- void GlobalFit();
- void ChiSquareSelection();
- void MergeTracks( vector<AliITSUTrackCooked> &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<Cell> 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
+++ /dev/null
-#ifndef ALIITSUTRACKERSAAUX_H
-#define ALIITSUTRACKERSAAUX_H
-
-#include <vector>
-using std::vector;
-#include <algorithm>
-using std::sort;
-#include "AliExternalTrackParam.h"
-#include "AliITSUTrackCooked.h"
-#include "AliITSUAux.h"
-
-#include <TMath.h>
-
-#include <Riostream.h>
-using std::cout;
-using std::endl;
-
-//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
-template <typename T>
-struct Comparison { //Adapted from TMath ROOT code
- Comparison(vector<T> *d) : fData(d) {}
-
- bool operator()(int i1, int i2) {
- return fData->at(i1).CmpVar() < fData->at(i2).CmpVar();
- }
- vector<T> *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<int> 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<SpacePoint> cmp(&Points);
- sort(Index.begin(),Index.end(),cmp);
- }
-
- void Clear() {
- Index.clear();
- Points.clear();
- N=0;
- }
-
- vector<int> Index;
- int N;
- vector<SpacePoint> 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 <AliITSUTrackCooked> {
- Comparison(vector<AliITSUTrackCooked> *d) : fData(d) {}
-
- bool operator()(int i1, int i2) {
- return (fData->at(i1).GetChi2() < fData->at(i2).GetChi2());
- }
- vector<AliITSUTrackCooked> *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
AliITSUClusterPix.cxx
AliITSUSeed.cxx
AliITSUTrackerGlo.cxx
- AliITSUTrackerSA.cxx
+ AliITSUCATracker.cxx
+ AliITSUCATrackingStation.cxx
AliITSUTrackerCooked.cxx
AliITSUTrackCooked.cxx
AliITSUTrackCond.cxx
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)
#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+;
#include <TGeoManager.h>
#include <TGeoGlobalMagField.h>
#include <TStopwatch.h>
+#include <TCanvas.h>
+#include <TVirtualPad.h>
#include "AliTracker.h"
#include "AliGeomManager.h"
#include "AliESDEvent.h"
#include "AliITSURecoParam.h"
#include "AliITSUReconstructor.h"
-#include "AliITSUTrackerSA.h"
+#include "AliITSUCATracker.h"
#endif
extern TSystem *gSystem;
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);
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; i<nEvents; i++) {
+ //for ( Int_t i=1; i<2; i++) {
cout<<"\nEvent number "<<i<<endl;
rl->GetEvent(i);
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();
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) {
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff