-//-------------------------------------------------------------------------
+//--------------------------------------------------------------------------------
// 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>
// Vc library
//#include "Vc/Vc"
-//#include "AliITSUTrackerSAauxVc.h" // Structs and other stuff using Vc library
+//#include "AliITSUTrackerSAauxVc.h" // Structs and other stuff using Vc library
#include "AliLog.h"
#include "AliESDEvent.h"
#include "AliITSUClusterPix.h"
#include <Riostream.h>
+using std::cout;
+using std::endl;
+using std::flush;
-
-//#include "AliITSUtrackSA.h" // Some dedicated SA track class ?
+//#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)
- //
- ,fClusters(),
- fClustersTC(),
- fDoublets(),
- fIndex(),
- fNClusters(),
- fNDoublets(),
- fPhiCut(0.05),
- fRPhiCut(0.03),
- fZCut(0.01)
+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
//--------------------------------------------------------------------
- for(Int_t i=0;i<7;++i) {
- fClusters[i].reserve(5000);
- }
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),
+AliITSUTrackerSA::AliITSUTrackerSA(const AliITSUTrackerSA &t): AliTracker(t),
+ fReconstructor(t.fReconstructor),
+ fITS(t.fITS),
+ fMatLUT(t.fMatLUT),
+ fUseMatLUT(t.fUseMatLUT),
+ fCurrMass(t.fCurrMass),
//
- fClusters(),
- fClustersTC(),
- fIndex(),
- fNClusters(),
- fNDoublets(),
- fPhiCut(),
- fRPhiCut(),
- fZCut()
+ fClustersTC(),
+ fChi2Cut(fgkChi2Cut),
+ fPhiCut(),
+ fRPhiCut(),
+ fZCut()
{
//--------------------------------------------------------------------
// The copy constructor is protected
}
//________________________________________________________________________________
-AliITSUTrackerSA::~AliITSUTrackerSA()
+AliITSUTrackerSA::~AliITSUTrackerSA()
{
// d-tor
delete fMatLUT;
// Possibly, increment the seeds with additional clusters (Kalman)
- // Possibly, (re)fit the found tracks
+ // Possibly, (re)fit the found tracks
- // Three iterations:
+ // Three iterations:
// - High momentum first;
- // - Low momentum with vertex constraint;
- // - Everything else.
+ // - Low momentum with vertex constraint;
+ // - Everything else.
- MakeDoublets(); // To be implemented
- //MakeTriplets(); // Are triplets really necessary? MFT does not use them.
- CASelection(event);
- //GlobalFit(); // To be implemented
- //ChiSquareSelection(); // To be implemented
+ CellsCreation(0);
+ CellularAutomaton(event);
+ // VertexFinding();
+ // CellsCreation(1);
+ // CellularAutomaton(event);
return 0;
}
//________________________________________________________________________________
-Int_t AliITSUTrackerSA::PropagateBack(AliESDEvent */*event*/) {
+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*/) {
+Int_t AliITSUTrackerSA::RefitInward(AliESDEvent * event) {
//--------------------------------------------------------------------
- // Some final refit, after the outliers get removed by the smoother ?
+ // 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,
//
for(int iL=0; iL<7; ++iL) {
fClustersTC[iL]=*fITS->GetLayerActive(iL)->GetClustersAddress();
- TClonesArray *clCont=fClustersTC[iL];
- fNClusters[iL]=clCont->GetEntriesFast();
- Float_t phi[fNClusters[iL]];
- fIndex[iL] = new Int_t[fNClusters[iL]];
- for(int iC=0;iC<fNClusters[iL];++iC) {
- const AliITSUClusterPix *cl = (AliITSUClusterPix*)clCont->At(iC);
+ 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);
- phi[iC] = pos[0]==0.f ? TMath::PiOver2() : TMath::ATan2(pos[1]-GetY(),pos[0]-GetX());
- float angle=0.f; // TO BE UNDERSTOOD: DO I STILL NEED THE STATION ANGLE IF I USE THE GLOBAL COVARIANCE MATRIX?
- fClusters[iL].push_back(itsCluster(pos[0],pos[1],pos[2],cl->GetSigmaY2(),cl->GetSigmaZ2(),cl->GetSigmaYZ(),phi[iC],angle));
+ 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);
+ }
}
- TMath::Sort(fNClusters[iL],phi,fIndex[iL],kFALSE);
+
+ fLayer[iL].Sort();
+
}
-// #ifdef __DEBUG__
-// PrintInfo("clusters");
-// #endif
return 0;
}
//--------------------------------------------------------------------
// This function unloads ITSU clusters from the RAM
//--------------------------------------------------------------------
- for(int i=0;i<7;++i) {
- fClusters[i].clear();
- fNClusters[i]=0;
- delete fIndex[i];
- }
- for(int i=0;i<6;++i) fDoublets[i].clear();
-}
+ for(int i=0;i<7;++i)
+ fLayer[i].Clear();
+ for(int i=0;i<5;++i)
+ fCells[i].clear();
-//________________________________________________________________________________
-AliCluster *AliITSUTrackerSA::GetCluster(Int_t /*index*/) const {
- //--------------------------------------------------------------------
- // Return pointer to a given cluster
- //--------------------------------------------------------------------
- return 0; // replace with an actual pointer
}
//________________________________________________________________________________
-void AliITSUTrackerSA::CASelection(AliESDEvent *event) {
+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
+ // 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.
-
- cout << "Begin of the CA selection" << endl;
- for( int iL = 1; iL < 6; ++iL ) {
-
- const itsCluster* clusters1 = &fClusters[iL-1][0];
- const itsCluster* clusters2 = &fClusters[iL][0];
- const itsCluster* clusters3 = &fClusters[iL+1][0];
-
- const nPlets* doublets1 = &fDoublets[iL-1][0];
- nPlets* doublets2 = &fDoublets[iL][0];
-
- for ( int iD2 = 0; iD2 < fNDoublets[iL]; ++iD2 ) {
- for ( int iD1 = 0; iD1 < fNDoublets[iL-1]; ++iD1 ) {
- const int id1 = doublets1[iD1].id0;
- const int id2 = doublets2[iD2].id0;
- if ( id1 == id2 ) {
- if ( doublets2[iD2].level <= ( doublets1[iD1].level + 1 ) ) {
- const int id3 = doublets2[iD2].id1;
- const float r3 = Sqrt( clusters3[id3].x * clusters3[id3].x + clusters3[id3].y * clusters3[id3].y );
- const float r2 = Sqrt( clusters2[id2].x * clusters2[id2].x + clusters2[id2].y * clusters2[id2].y );
- const float extrZ3 = doublets1[iD1].tanLambda * ( r3 - r2 ) + clusters2[id2].z ;
- //cout << extrZ3 << " " << clusters3[id3].z << " " << Abs ( extrZ3 - clusters3[id3].z ) << endl;
- if ( Abs ( extrZ3 - clusters3[id3].z ) < fZCut ) {
- //cout << "OK Z doublets: "<< iL-1 << "," << iD1 << "\t" << iL << "," <<iD2 << endl;
- const float det = (clusters1[id1].x - GetX())*(clusters2[id2].y - GetY()) - (clusters2[id2].x-GetX() )*(clusters1[id1].y - GetY()); // (GetX() - clusters2[id2].x)*(clusters1[id1].y - clusters2[id2].y) - (GetY() - clusters2[id2].y)*(clusters1[id1].x - clusters2[id2].x);
- //cout << det << endl;
- if ( Abs(det) <= 1e-12 ) {
- // linear extrapolation to next layer
- const float dsq = ( doublets1[iD1].tanPhi * (clusters3[id3].x + clusters2[id2].x) + clusters3[id3].y - clusters2[id2].y ) *
- ( doublets1[iD1].tanPhi * (clusters3[id3].x + clusters2[id2].x) + clusters3[id3].y - clusters2[id2].y ) / (1 + doublets1[iD1].tanPhi * doublets1[iD1].tanPhi );
- if ( dsq < fRPhiCut*fRPhiCut ) {
- doublets2[iD2].level += doublets1[iD1].level;
- doublets2[iD2].neighbours.push_back(iD1);
- }
- } else {
- const float r1sq = clusters1[id1].x * clusters1[id1].x + clusters1[id1].y * clusters1[id1].y ;
- const float rvsq = GetX() * GetX() + GetY() * GetY();
- const float deta = (rvsq - r1sq) * (clusters2[id2].y - GetY()) - (rvsq - r2*r2) * (clusters1[id1].y - GetY());
- const float detb = - (rvsq - r1sq) * (clusters2[id2].x - GetX()) + (rvsq - r2*r2) * (clusters1[id1].x - GetX()) ;
- const float a = deta/det ;
- const float b = detb/det ;
- const float c = -rvsq - a * GetX() - b * GetY();
- const float rc = Sqrt( a*a/4.f + b*b/4.f - c );
- const float d = Sqrt( (a/2.f + clusters3[id3].x) * (a/2.f + clusters3[id3].x) + (b/2.f + clusters3[id3].y) * (b/2.f + clusters3[id3].y) );
- //cout << d << " " << rc << " " << d - rc << endl;
- if ( Abs( d - rc ) < fRPhiCut ) {
- doublets2[iD2].level += doublets1[iD1].level;
- doublets2[iD2].neighbours.push_back(iD1);
- }
- }
- }
- }
- }
+ 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);
+ }
}
}
}
-
- /*#ifdef __DEBUG__
- PrintInfo("doublets");
- #endif*/
-
- // Hic sunt leones: the following code could be optimised to be iterative. But now I don't have time.
- for ( int level = 6; level >= 2 ; --level ) {
- vector<trackC> candidates;
- int nCand=0;
- const int limit = level > 5 ? 5 : level;
- for ( int doubl = 5; doubl >= limit; --doubl ) {
- for ( int iD = 0; iD < fNDoublets[doubl]; ++iD ) {
- if ( fDoublets[doubl][iD].level == level ) {
- cout << "Pushing new candidate" << endl;
- candidates.push_back(trackC());
- candidates[nCand].fPoints[doubl*2+2] = fDoublets[doubl][iD].id1;
- candidates[nCand].fPoints[doubl*2] = fDoublets[doubl][iD].id0;
- CandidatesTreeTraversal(candidates, iD, doubl, nCand);
- }
+
+ 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;
}
}
- nCand = candidates.size();
- Double_t rDest = fITS->GetRMax();
- cout << "Number of candidates at level " << level << ": " << candidates.size() << endl;
- int index[nCand];
- for( int i = 0; i < nCand; ++i ) index[i] = i;
- for ( int cand = 0; cand < nCand; ++cand ) {
- int count = 0;
- for ( int j = 0; j < 14; ++j ) {
- if( candidates[cand].fPoints[j] > -1 ) ++count;
+ // 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;
}
- cout << "Candidates " << cand << ", number of points: " << count << endl;
- // candidates[cand].ResetCovariance(50.);
- InitTrackParams(candidates[cand]);
- cout << "Fit cnd: " << cand << " " << candidates[cand] << endl;
- candidates[cand].fChi2 = RefitTrack( (AliExternalTrackParam*)&candidates[cand], candidates[cand].fPoints, rDest ,-1);
+
+ 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);
}
- CompDesc comp(&candidates);
- sort(index,index+nCand,comp);
- for ( int cand = 0; cand < nCand; ++cand ) {
+
+ 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];
-#ifdef __DEBUG__
- //cout << ii << " " << candidates[ii] << endl;
-#endif
- if ( candidates[ii].fChi2 < 0. ) break;
+
+ 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 ( unsigned int point = 0; point < 14; ++point ) {
- if ( candidates[ii].fPoints[point] != -1 ) {
- if( !(fClusters[ point/2 ][ candidates[ii].fPoints[point] ].isUsed ) ) {
- fClusters[ point/2 ][ candidates[ii].fPoints[point] ].isUsed = true;
- } else {
- goodTrack = false;
- }
- }
+ 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 ) {
- AliESDtrack outTrack;
- outTrack.SetOuterParam((AliExternalTrackParam*)&candidates[ii],AliESDtrack::kITSpureSA);
- event->AddTrack(&outTrack);
+ 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++;
- // A first try...
- // for ( int level = 6; level >=3 ; --level ) {
- // vector<int> points[7];
- // for ( int doubl = 5; doubl >= 2; --doubl ) {
- // if( ( doubl + 1 - level ) < 0 ) break;
- // for ( int iD = 0; iD < fNDoublets[doubl]; ++iD ) {
- // if ( fDoublets[doubl][iD].level == level ) {
- // points[doubl+1].push_back( fDoublets[doubl][iD].id1 );
- // points[doubl].push_back( fDoublets[doubl][iD].id0 );
-
- // vector<int> currentVector = (fDoublets[doubl][iD].neighbours);
- // for( int iL = 1; iL <= doubl; ++iL ) {
- // const nPlets * doublets = &fDoublets[doubl-iL][0];
- // vector<int> tmp;
- // for ( unsigned int iV = 0 ; iV < currentVector.size() ; ++iV ) {
- // points[doubl-iL].push_back( doublets[ currentVector[ iV ] ].id0 ) ;
- // for ( unsigned int iN = 0 ; iN < doublets[ currentVector[ iV ] ].neighbours.size(); ++iN ) {
- // tmp.push_back( doublets[ currentVector[ iV ] ].neighbours[ iN ] );
- // }
- // }
- // currentVector.swap( tmp );
- // }
- // }
- // }
- // }
- // }
+ // 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::MakeDoublets() {
- // Make associations between two points on adjacent layers within an azimuthal window.
+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
- cout << "Vertex of used by the tracker: " << GetX() << " " << GetY() << " " << GetZ() << endl;
+ float phiCut = 7.f;
+ if( cutLevel==0 ) phiCut = fPhiCut;
+ // Doublets creation
+ vector<Cell> doublets[6];
for( int iL = 0 ; iL < 6 ; ++iL ) {
- fNDoublets[iL] = 0;
- const itsCluster* clusters1 = &fClusters[iL][0];
- const itsCluster* clusters2 = &fClusters[iL+1][0];
-
- // 0 - 2Pi junction treatment (part I)
- for ( int iCC1 = 0 ; iCC1 < fNClusters[iL] ; ++iCC1 ) {
- bool flag = true;
- const int iC1 = fIndex[iL][iCC1];
- for ( int iCC2 = fNClusters[iL+1]-1; iCC2 > -1 ; --iCC2 ) {
- const int iC2 = fIndex[iL+1][iCC2];
- if( (TMath::TwoPi() - (clusters2[iC2].phi-clusters1[iC1].phi) ) < fPhiCut ) {
- fDoublets[iL].push_back(nPlets(iC1,iC2));
- fDoublets[iL][fNDoublets[iL]].tanPhi = (clusters1[iC1].y-clusters2[iC2].y)/(clusters1[iC1].x-clusters2[iC2].x);
- float r1 = Sqrt(clusters1[iC1].x * clusters1[iC1].x + clusters1[iC1].y * clusters1[iC1].y);
- cout << clusters2[iC2].x * clusters2[iC2].x + clusters2[iC2].y * clusters2[iC2].y << flush << endl;
- float r2 = Sqrt(clusters2[iC2].x * clusters2[iC2].x + clusters2[iC2].y * clusters2[iC2].y);
- fDoublets[iL][fNDoublets[iL]].tanLambda = (clusters1[iC1].z-clusters2[iC2].z)/(r1-r2);
- ++fNDoublets[iL];
- flag = false;
- } else break;
-
- }
- if (flag) break;
- }
-
-
- // "Central" points
- for ( int iCC1 = 0 ; iCC1 < fNClusters[iL] ; ++iCC1 ) {
- const int iC1 = fIndex[iL][iCC1];
- for ( int iCC2 = 0; iCC2 < fNClusters[iL+1] ; ++iCC2 ) {
- const int iC2 = fIndex[iL+1][iCC2];
- if( Abs( clusters1[iC1].phi - clusters2[iC2].phi ) < fPhiCut ) {
- fDoublets[iL].push_back(nPlets(iC1,iC2));
- fDoublets[iL][fNDoublets[iL]].tanPhi = (clusters1[iC1].y-clusters2[iC2].y)/(clusters1[iC1].x-clusters2[iC2].x);
- float r1 = Sqrt(clusters1[iC1].x * clusters1[iC1].x + clusters1[iC1].y * clusters1[iC1].y);
- float r2 = Sqrt(clusters2[iC2].x * clusters2[iC2].x + clusters2[iC2].y * clusters2[iC2].y);
- fDoublets[iL][fNDoublets[iL]].tanLambda = (clusters1[iC1].z-clusters2[iC2].z)/(r1-r2);
- ++fNDoublets[iL];
- } else if( clusters2[iC2].phi - clusters1[iC1].phi > fPhiCut ) break;
-
+ 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;
}
}
+ }
- // 0 - 2Pi junction treatment (part II)
- for ( int iCC1 = fNClusters[iL]-1; iCC1 > -1 ; --iCC1 ) {
- bool flag = true;
- const int iC1 = fIndex[iL][iCC1];
- for ( int iCC2 = 0; iCC2 < fNClusters[iL+1] ; ++iCC2 ) {
- const int iC2 = fIndex[iL+1][iCC2];
- if( (TMath::TwoPi() - (clusters1[iC1].phi-clusters2[iC2].phi) ) < fPhiCut ) {
- fDoublets[iL].push_back(nPlets(iC1,iC2));
- fDoublets[iL][fNDoublets[iL]].tanPhi = (clusters1[iC1].y-clusters2[iC2].y)/(clusters1[iC1].x-clusters2[iC2].x);
- float r1 = Sqrt(clusters1[iC1].x * clusters1[iC1].x + clusters1[iC1].y * clusters1[iC1].y);
- float r2 = Sqrt(clusters2[iC2].x * clusters2[iC2].x + clusters2[iC2].y * clusters2[iC2].y);
- fDoublets[iL][fNDoublets[iL]].tanLambda = (clusters1[iC1].z-clusters2[iC2].z)/(r1-r2);
- ++fNDoublets[iL];
- flag = false;
- } else 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;
+ }
+ }
}
-
- if (flag) break;
}
-
}
- // #ifdef __DEBUG__
- // PrintInfo("doublets");
- // #endif
+
}
+
//______________________________________________________________________________
-Bool_t AliITSUTrackerSA::InitTrackParams(trackC &track)
+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();i--;) if (track.fPoints[i<<0x1]>-1) lrOcc[nCl++] = i;
+ 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));
- cout << track << endl;
return kFALSE;
}
//
- int lr0 = lrOcc[0];
- int lr1 = lrOcc[nCl/2];
- int lr2 = lrOcc[nCl-1];
+ const int lr0 = lrOcc[0];
+ const int lr1 = lrOcc[nCl/2];
+ const int lr2 = lrOcc[nCl-1];
//
- const itsCluster& cl0 = fClusters[lr0][ track.fPoints[lr0<<0x1] ];
- const itsCluster& cl1 = fClusters[lr1][ track.fPoints[lr1<<0x1] ];
- const itsCluster& cl2 = fClusters[lr2][ track.fPoints[lr2<<0x1] ];
- double cv = Curvature(cl0.x,cl0.y, cl1.x,cl1.y, cl2.x,cl2.y);
- double tgl = (cl2.z-cl0.z)/TMath::Sqrt((cl2.x-cl0.x)*(cl2.x-cl0.x)+(cl2.y-cl0.y)*(cl2.y-cl0.y));
- // double phi = TMath::ATan2((cl2.y-cl1.y),(cl2.x-cl1.x));
+ 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( track.fPoints[lr0<<0x1] );
+ 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 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
+ 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<trackC> &track, int &iD, int &doubl, int &nCand) {
-
- if ( doubl < 1 ) {
-#ifdef __DEBUG__
- cout << "ERROR IN CandidatesTreeTraversal" << endl;
-#endif
- return;
- }
-
- int currD = doubl-1;
- track[nCand].fPoints[ 2 * currD ] = fDoublets[ currD ] [ fDoublets[ doubl ][ iD ].neighbours[0] ].id0;
+void AliITSUTrackerSA::CandidatesTreeTraversal(vector<Road> &candidates, const int &iD, const int &doubl) {
+ if ( doubl < 0 ) return;
- if ( fDoublets[ currD ] [ fDoublets[ doubl ][ iD ].neighbours[0] ].level == 1 ) {
-
- cout << "Setting track " << nCand << " information using its first cluster: " << endl;
- cout << "Layer " << currD << ", id " << track[nCand].fPoints[ 2 * currD ] <<endl;
- AliITSUClusterPix* clus = (AliITSUClusterPix*)fClustersTC[ currD ]->At( track[nCand].fPoints[ 2 * currD ] ); // assign your 1st cluster (from which the fit is starting)
- AliITSURecoLayer* lr = fITS->GetLayerActive(currD) ; // assign the layer which the cluster belongs to
+ 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];
- AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
- double x = sens->GetXTF() + clus->GetX();
- double alp = sens->GetPhiTF();
- double par[5] = {clus->GetY(),clus->GetZ(),0,0,0.01};
- 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[nCand].Set(x,alp,par,cov);
- // cout << "after set " << track[nCand] << endl;
- return;
- } else {
- CandidatesTreeTraversal(track,fDoublets[doubl][iD].neighbours[0],currD,nCand);
- }
-
- for( unsigned int iD2 = 1 ; iD2 < fDoublets[ doubl ][ iD ].neighbours.size() ; ++iD2 ) {
- track.push_back(track[nCand++]);
- track[nCand].fPoints[ 2 * currD ] = fDoublets[ currD ] [ fDoublets[ doubl ][ iD ].neighbours[ iD2 ] ].id0;
- if ( fDoublets[ currD ] [ fDoublets[ doubl ][ iD ].neighbours[iD2] ].level == 1 ) {
- AliITSUClusterPix* clus = (AliITSUClusterPix*)fClustersTC[ currD ]->At( track[nCand].fPoints[ 2 * currD ] ); // assign your 1st cluster (from which the fit is starting)
- AliITSURecoLayer* lr = fITS->GetLayerActive(currD) ; // assign the layer which the cluster belongs to
-
- AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
- double x = sens->GetXTF() + clus->GetX();
- double alp = sens->GetPhiTF();
- double par[5] = {clus->GetY(),clus->GetZ(),0,0,0};
- 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[nCand].Set(x,alp,par,cov);
- return;
- } else {
- CandidatesTreeTraversal(track,fDoublets[doubl][iD].neighbours[iD2],currD,nCand);
+ 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)
+ 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
dir = rCurr<rDest ? 1 : -1;
lrStart = fITS->FindFirstLayerID(rCurr,dir);
lrStop = fITS->FindLastLayerID(rDest,dir); // lr id before which we have to stop
- cout << "Start/End layers and direction " << lrStart << " " << lrStop << " " << dir << endl;
//
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));
+ "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++;}
int lrStop1 = lrStop+dir;
for (int ilr=lrStart;ilr!=lrStop1;ilr+=dir) {
AliITSURecoLayer* lr = fITS->GetLayer(ilr);
- if ( dir*(rCurr-lr->GetR(dir))>0) {
- cout << ilr << " passed!" << endl;
- continue;
- } // this layer is already passed
+ 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) {
- cout << ilr << " is inactive or without cluster for current candidates" << endl;
- continue;
- }
- cout << "OK layer " << ilr << endl;
+ 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;
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())) { cout << "failed rotation" << endl; return -1; }
+ if (!tmpTr.Rotate(sens->GetPhiTF())) return -1;
//
double xClus = sens->GetXTF()+clus->GetX();
if (!transportedToLayer) {
- if (ilr!=lrStart && !TransportToLayerX(&tmpTr,lrStart,ilr,xClus)) {
- cout << "failed transport to the entrance" << endl; return -1; } // go to the entrance to the layer
- lrStart = ilr;
- transportedToLayer = kTRUE;
+ 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)) { cout << "failed propagation of the seed X:" << xClus << endl; tmpTr.Print(); return -1; }
+ 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) ) { cout << "failed update of the covariance" << endl; return -1; }
+ if ( !tmpTr.Update(p,cov) ) return -1;
if (++nclFit==nCl && stopCond<0) {
- *trc = tmpTr;
- // printf("Fit chi2: %f for %d clusters\n",chi2,nclFit);
- return chi2; // it was requested to not propagate after last update
+ *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.
//
- // printf("Fit chi2: %f for %d clusters\n",chi2,nclFit);
- //
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
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;
+ if (!seed->PropagateTo(x,bz)) return kFALSE;
double ds = 0;
if (matCorr || updTime) {
xyz0[0]=xyz1[0]; // global pos at the beginning of step
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;
+ 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);
+ 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);
double xToGo = lrTo->GetR(-dir);
if (rLim>0) {
if (dir>0) {
- if (rLim<xToGo) {xToGo = rLim; limReached = kTRUE;}
+ if (rLim<xToGo) {xToGo = rLim; limReached = kTRUE;}
}
else {
- if (rLim>xToGo) {xToGo = rLim; limReached = kTRUE;}
+ if (rLim>xToGo) {xToGo = rLim; limReached = kTRUE;}
}
}
// double xts = xToGo;
}
//____________________________________________________________________
-Double_t AliITSUTrackerSA::Curvature(Double_t x1,Double_t y1,Double_t
-x2,Double_t y2,Double_t x3,Double_t y3)
-{
+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
+ //calculates the curvature of track
Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
- if(den==0) return 0;
+ 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.;
+ Double_t xc= -a/2.;
- if((a*a+b*b-4*c)<0) return 0;
+ if((a*a+b*b-4*c)<0) return 0.;
Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
- if(rad==0) return 0;
-
+ 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;
-
+ // if((x1<0 && y1<0 && x1<xc)) rad*=-1;
+ // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
+
return 1/rad;
-
-}
-#ifdef __DEBUG__
-//____________________________________________________
-void AliITSUTrackerSA::PrintInfo(TString what) {
- //
- if( what.Contains("clusters") ) {
- cout << "Dumping clusters info" << endl;
- for ( int i = 0; i < 7; ++i ) {
- cout << "**** Layer " << i << " ****" << endl;
- for ( int c = 0; c < fNClusters[i]; ++c ) {
- cout << "*** Cluster " << c << " ***" <<endl;
- cout << fClusters[i][c] << endl;
- }
- }
- }
- //
- if( what.Contains("doublets") ) {
- cout << "Dumping doublets info" << endl;
- for ( int i = 0; i < 6; ++i ) {
- cout << "**** Doublets array " << i << " ****" << endl;
- for ( int c = 0; c < fNDoublets[i]; ++c ) {
- cout << "*** Doublet " << c << " ***" <<endl;
- cout << fDoublets[i][c] << endl;
- }
- }
- }
}
-#endif
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff