[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUTrackerSA.cxx

//-------------------------------------------------------------------------
//               Implementation of the ITS tracker class
//    It reads AliITSUClusterPix clusters and and fills the ESD with tracks
//-------------------------------------------------------------------------

#include <TBranch.h>
#include <TMath.h>
using TMath::Abs;
using TMath::Sort;
using TMath::Sqrt;
#include <TTree.h>

// Vc library
//#include "Vc/Vc"
//#include "AliITSUTrackerSAauxVc.h" // Structs and other stuff using Vc library  

#include "AliESDEvent.h"
#include "AliITSUClusterPix.h"
#include "AliITSUTrackerSA.h"

//#include "AliITSUtrackSA.h"      // Some dedicated SA track class ?  

ClassImp(AliITSUTrackerSA)

//________________________________________________________________________________
AliITSUTrackerSA::AliITSUTrackerSA() : AliTracker(), 
  fClusters(),
  fClustersTC(),
  fDoublets(),
  fIndex(),
  fNClusters(),
  fNDoublets(),
  fPhiCut(0.05),
  fRPhiCut(0.01),
  fZCut(0.005)
{
  //--------------------------------------------------------------------
  // This default constructor needs to be provided
  //--------------------------------------------------------------------
  for(Int_t i=0;i<7;++i) {
    fClusters[i].reserve(5000);
  }
}

//________________________________________________________________________________
AliITSUTrackerSA::AliITSUTrackerSA(const AliITSUTrackerSA &t): 
  AliTracker(t),
  fClusters(),
  fClustersTC(),
  fIndex(),
  fNClusters(),
  fNDoublets(),
  fPhiCut(),
  fRPhiCut(),
  fZCut()
{
  //--------------------------------------------------------------------
  // The copy constructor is protected
  //--------------------------------------------------------------------
}

//________________________________________________________________________________
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 

  // Tree iterations: 
  // - High momentum first;
  // - Low momentum with vertex constraint; 
  // - Everything else. 

  MakeDoublets();       // To be implemented
  //MakeTriplets();       // Are triplets really necessary? MFT does not use them. 
  CASelection();        // To be implemented
  GlobalFit();          // To be implemented
  ChiSquareSelection(); // To be implemented

  return 0;
}

//________________________________________________________________________________
Int_t AliITSUTrackerSA::PropagateBack(AliESDEvent */*event*/) {
  //--------------------------------------------------------------------
  // Here, we implement the Kalman smoother ?
  // The clusters must be already loaded
  //--------------------------------------------------------------------

  return 0;
}

//________________________________________________________________________________
Int_t AliITSUTrackerSA::RefitInward(AliESDEvent */*event*/) {
  //--------------------------------------------------------------------
  // Some final refit, after the outliers get removed by the smoother ?  
  // The clusters must be loaded
  //--------------------------------------------------------------------

  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
  //--------------------------------------------------------------------
  
  for(Int_t i=0;i<7;++i) {
    TBranch* br = cluTree->GetBranch(Form("ITSRecPoints%d",i));
    if (!br) return -1;
    br->SetAddress(&fClustersTC[i]);
  }    
  cluTree->GetEntry(0);
  
  for(int iL=0; iL<7; ++iL) {
    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);
      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));
    }
    TMath::Sort(fNClusters[iL],phi,fIndex[iL],kFALSE);
  }
  
  return 0;
}

//________________________________________________________________________________
void AliITSUTrackerSA::UnloadClusters() {
  //--------------------------------------------------------------------
  // 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();
}

//________________________________________________________________________________
AliCluster *AliITSUTrackerSA::GetCluster(Int_t /*index*/) const {
  //--------------------------------------------------------------------
  //       Return pointer to a given cluster
  //--------------------------------------------------------------------
  return 0;  // replace with an actual pointer 
}

//________________________________________________________________________________
void AliITSUTrackerSA::CASelection() {
  // 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.

  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].id1;
	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( clusters3[id2].x * clusters2[id2].x + clusters2[id2].y * clusters2[id2].y );
	    const float extrZ3 = doublets1[iD1].tanLambda * ( r3 - r2 ) + clusters3[id3].z ;
	    if ( Abs ( extrZ3 - clusters3[id3].z ) < fZCut ) {
	      const float det = (GetX() - clusters2[id2].x)*(GetY() - clusters1[id1].y) - (GetY() - clusters2[id2].y)*(GetX() - clusters1[id1].x);
	      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 = (r2*r2-rvsq) * (GetY() - clusters1[id1].y) - (r1sq-rvsq) * (GetY() - clusters2[id2].y);
		const float detb =  (r2*r2-rvsq) * (GetX() - clusters1[id1].x) - (r1sq-rvsq) * (GetX() - clusters2[id2].x) ;
		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) );
		if ( ( d - rc ) < fRPhiCut ) {
		  doublets2[iD2].level += doublets1[iD1].level;
		  doublets2[iD2].neighbours.push_back(iD1);
		}
	      }
	    }
	  }
	}
      }
    }
  }

  for ( int level = 6; level >=3 ; --level ) {
    //vector<int> points[6];
    for ( int doubl = 5; doubl >= 0; --doubl ) {
      if( ( doubl + 1 - level ) < 0 ) break;
      for ( int iD = 0; iD < fNDoublets[doubl]; ++iD ) {
	if ( fDoublets[doubl][iD].level == level ) {
	  
	}
      }
    }
  }
}

//________________________________________________________________________________
void AliITSUTrackerSA::MakeDoublets() {
  // 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

  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 iC1 = 0 ; iC1 < fNClusters[iL] ; ++iC1 ) {
      bool flag = true;
      for ( int iC2 = fNClusters[iL]-1; iC2 > -1 ; --iC2 ) {
	
	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);
	  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 iC1 = 0 ; iC1 < fNClusters[iL] ; ++iC1 ) {

      for ( int iC2 = 0; iC2 < fNClusters[iL+1] ; ++iC2 ) {
	
	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;
      
      }

    }

    // 0 - 2Pi junction treatment (part II)
    for ( int iC1 = fNClusters[iL]-1; iC1 > -1 ; --iC1 ) {
      bool flag = true;

      for ( int iC2 = 0; iC2 < fNClusters[iL+1] ; ++iC2 ) {

	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;

      }

      if (flag) break;
    }
  
  }

}
Commit	Line	Data
2cbb5f31	1	//-------------------------------------------------------------------------
	2	// Implementation of the ITS tracker class
	3	// It reads AliITSUClusterPix clusters and and fills the ESD with tracks
	4	//-------------------------------------------------------------------------
	5
	6	#include <TBranch.h>
	7	#include <TMath.h>
	8	using TMath::Abs;
	9	using TMath::Sort;
	10	using TMath::Sqrt;
	11	#include <TTree.h>
	12
	13	// Vc library
	14	//#include "Vc/Vc"
	15	//#include "AliITSUTrackerSAauxVc.h" // Structs and other stuff using Vc library
	16
	17	#include "AliESDEvent.h"
	18	#include "AliITSUClusterPix.h"
	19	#include "AliITSUTrackerSA.h"
	20
	21	//#include "AliITSUtrackSA.h" // Some dedicated SA track class ?
	22
	23	ClassImp(AliITSUTrackerSA)
	24
	25	//________________________________________________________________________________
	26	AliITSUTrackerSA::AliITSUTrackerSA() : AliTracker(),
	27	fClusters(),
	28	fClustersTC(),
	29	fDoublets(),
	30	fIndex(),
	31	fNClusters(),
	32	fNDoublets(),
	33	fPhiCut(0.05),
	34	fRPhiCut(0.01),
	35	fZCut(0.005)
	36	{
	37	//--------------------------------------------------------------------
	38	// This default constructor needs to be provided
	39	//--------------------------------------------------------------------
	40	for(Int_t i=0;i<7;++i) {
	41	fClusters[i].reserve(5000);
	42	}
	43	}
	44
	45	//________________________________________________________________________________
	46	AliITSUTrackerSA::AliITSUTrackerSA(const AliITSUTrackerSA &t):
	47	AliTracker(t),
	48	fClusters(),
	49	fClustersTC(),
	50	fIndex(),
	51	fNClusters(),
	52	fNDoublets(),
	53	fPhiCut(),
	54	fRPhiCut(),
	55	fZCut()
	56	{
	57	//--------------------------------------------------------------------
	58	// The copy constructor is protected
	59	//--------------------------------------------------------------------
	60	}
	61
	62	//________________________________________________________________________________
	63	Int_t AliITSUTrackerSA::Clusters2Tracks(AliESDEvent /event*/) {
	64	//--------------------------------------------------------------------
65	// This is the main tracking function
66	// The clusters must already be loaded
67	//--------------------------------------------------------------------
68
69	// Possibly, create the track "seeds" (combinatorial)
70
71	// Possibly, increment the seeds with additional clusters (Kalman)
72
73	// Possibly, (re)fit the found tracks
74
75	// Tree iterations:
76	// - High momentum first;
77	// - Low momentum with vertex constraint;
78	// - Everything else.
79
80	MakeDoublets(); // To be implemented
81	//MakeTriplets(); // Are triplets really necessary? MFT does not use them.
82	CASelection(); // To be implemented
83	GlobalFit(); // To be implemented
84	ChiSquareSelection(); // To be implemented
85
86	return 0;
87	}
88
89	//________________________________________________________________________________
90	Int_t AliITSUTrackerSA::PropagateBack(AliESDEvent /event*/) {
91	//--------------------------------------------------------------------
92	// Here, we implement the Kalman smoother ?
93	// The clusters must be already loaded
94	//--------------------------------------------------------------------
95
96	return 0;
97	}
98
99	//________________________________________________________________________________
100	Int_t AliITSUTrackerSA::RefitInward(AliESDEvent /event*/) {
101	//--------------------------------------------------------------------
102	// Some final refit, after the outliers get removed by the smoother ?
103	// The clusters must be loaded
104	//--------------------------------------------------------------------
105
106	return 0;
107	}
108
109	Int_t AliITSUTrackerSA::LoadClusters(TTree *cluTree) {
110	//--------------------------------------------------------------------
111	// This function reads the ITSU clusters from the tree,
112	// sort them, distribute over the internal tracker arrays, etc
113	//--------------------------------------------------------------------
114
115	for(Int_t i=0;i<7;++i) {
116	TBranch* br = cluTree->GetBranch(Form("ITSRecPoints%d",i));
117	if (!br) return -1;
118	br->SetAddress(&fClustersTC[i]);
119	}
120	cluTree->GetEntry(0);
121
122	for(int iL=0; iL<7; ++iL) {
123	TClonesArray *clCont=&fClustersTC[iL];
124	fNClusters[iL]=clCont->GetEntriesFast();
125	Float_t phi[fNClusters[iL]];
126	fIndex[iL] = new Int_t[fNClusters[iL]];
127	for(int iC=0;iC<fNClusters[iL];++iC) {
128	const AliITSUClusterPix cl = (AliITSUClusterPix)clCont->At(iC);
129	float pos[3];
130	cl->GetGlobalXYZ(pos);
131	phi[iC] = pos[0]==0.f ? TMath::PiOver2() : TMath::ATan2(pos[1]-GetY(),pos[0]-GetX());
132	float angle=0.f; // TO BE UNDERSTOOD: DO I STILL NEED THE STATION ANGLE IF I USE THE GLOBAL COVARIANCE MATRIX?
133	fClusters[iL].push_back(itsCluster(pos[0],pos[1],pos[2],cl->GetSigmaY2(),cl->GetSigmaZ2(),cl->GetSigmaYZ(),phi[iC],angle));
134	}
135	TMath::Sort(fNClusters[iL],phi,fIndex[iL],kFALSE);
136	}
137
138	return 0;
139	}
140
141	//________________________________________________________________________________
142	void AliITSUTrackerSA::UnloadClusters() {
143	//--------------------------------------------------------------------
144	// This function unloads ITSU clusters from the RAM
145	//--------------------------------------------------------------------
146	for(int i=0;i<7;++i) {
147	fClusters[i].clear();
148	fNClusters[i]=0;
149	delete fIndex[i];
150	}
151	for(int i=0;i<6;++i) fDoublets[i].clear();
152	}
153
154	//________________________________________________________________________________
155	AliCluster AliITSUTrackerSA::GetCluster(Int_t /index*/) const {
156	//--------------------------------------------------------------------
157	// Return pointer to a given cluster
158	//--------------------------------------------------------------------
159	return 0; // replace with an actual pointer
160	}
161
162	//________________________________________________________________________________
163	void AliITSUTrackerSA::CASelection() {
164	// Here it's implemented the Cellular Automaton routine
165	// Firstly the level of each doublet is set according to the level of
166	// the neighbour doublets.
167	// Doublet are considered to be neighbour if they share one point and the
168	// phi and theta direction difference of the two is below a cut value.
169
170	for( int iL = 1; iL < 6; ++iL ) {
171
172	const itsCluster* clusters1 = &fClusters[iL-1][0];
173	const itsCluster* clusters2 = &fClusters[iL][0];
174	const itsCluster* clusters3 = &fClusters[iL+1][0];
175
176	const nPlets* doublets1 = &fDoublets[iL-1][0];
177	nPlets* doublets2 = &fDoublets[iL][0];
178
179	for ( int iD2 = 0; iD2 < fNDoublets[iL]; ++iD2 ) {
180	for ( int iD1 = 0; iD1 < fNDoublets[iL-1]; ++iD1 ) {
181	const int id1 = doublets1[iD1].id1;
182	const int id2 = doublets2[iD2].id0;
183	if ( id1 == id2 ) {
184	if ( doublets2[iD2].level <= ( doublets1[iD1].level + 1 ) ) {
185	const int id3 = doublets2[iD2].id1;
186	const float r3 = Sqrt( clusters3[id3].x * clusters3[id3].x + clusters3[id3].y * clusters3[id3].y );
187	const float r2 = Sqrt( clusters3[id2].x * clusters2[id2].x + clusters2[id2].y * clusters2[id2].y );
188	const float extrZ3 = doublets1[iD1].tanLambda * ( r3 - r2 ) + clusters3[id3].z ;
189	if ( Abs ( extrZ3 - clusters3[id3].z ) < fZCut ) {
190	const float det = (GetX() - clusters2[id2].x)(GetY() - clusters1[id1].y) - (GetY() - clusters2[id2].y)(GetX() - clusters1[id1].x);
191	if ( Abs(det) <= 1e-12 ) {
192	// linear extrapolation to next layer
193	const float dsq = ( doublets1[iD1].tanPhi * (clusters3[id3].x + clusters2[id2].x) + clusters3[id3].y - clusters2[id2].y ) *
194	( doublets1[iD1].tanPhi * (clusters3[id3].x + clusters2[id2].x) + clusters3[id3].y - clusters2[id2].y ) / (1 + doublets1[iD1].tanPhi * doublets1[iD1].tanPhi );
195	if ( dsq < fRPhiCut*fRPhiCut ) {
196	doublets2[iD2].level += doublets1[iD1].level;
197	doublets2[iD2].neighbours.push_back(iD1);
198	}
199	} else {
200	const float r1sq = clusters1[id1].x * clusters1[id1].x + clusters1[id1].y * clusters1[id1].y ;
201	const float rvsq = GetX() * GetX() + GetY() * GetY();
202	const float deta = (r2r2-rvsq) (GetY() - clusters1[id1].y) - (r1sq-rvsq) * (GetY() - clusters2[id2].y);
203	const float detb = (r2r2-rvsq) (GetX() - clusters1[id1].x) - (r1sq-rvsq) * (GetX() - clusters2[id2].x) ;
204	const float a = deta/det ;
205	const float b = detb/det ;
206	const float c = -rvsq - a * GetX() - b * GetY();
207	const float rc = Sqrt( aa/4.f + bb/4.f - c );
208	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) );
209	if ( ( d - rc ) < fRPhiCut ) {
210	doublets2[iD2].level += doublets1[iD1].level;
211	doublets2[iD2].neighbours.push_back(iD1);
212	}
213	}
214	}
215	}
216	}
217	}
218	}
219	}
220
221	for ( int level = 6; level >=3 ; --level ) {
222	//vector<int> points[6];
223	for ( int doubl = 5; doubl >= 0; --doubl ) {
224	if( ( doubl + 1 - level ) < 0 ) break;
225	for ( int iD = 0; iD < fNDoublets[doubl]; ++iD ) {
226	if ( fDoublets[doubl][iD].level == level ) {
227
228	}
229	}
230	}
231	}
232	}
233
234	//________________________________________________________________________________
235	void AliITSUTrackerSA::MakeDoublets() {
236	// Make associations between two points on adjacent layers within an azimuthal window.
237	// Under consideration:
238	// - track parameter estimation using the primary vertex position
239	// To do:
240	// - last iteration
241
242	for( int iL = 0 ; iL < 6 ; ++iL ) {
243	fNDoublets[iL] = 0;
244	const itsCluster* clusters1 = &fClusters[iL][0];
245	const itsCluster* clusters2 = &fClusters[iL+1][0];
246
247	// 0 - 2Pi junction treatment (part I)
248	for ( int iC1 = 0 ; iC1 < fNClusters[iL] ; ++iC1 ) {
249	bool flag = true;
250	for ( int iC2 = fNClusters[iL]-1; iC2 > -1 ; --iC2 ) {
251
252	if( (TMath::TwoPi() - (clusters2[iC2].phi-clusters1[iC1].phi) ) < fPhiCut ) {
253	fDoublets[iL].push_back(nPlets(iC1,iC2));
254	fDoublets[iL][fNDoublets[iL]].tanPhi = (clusters1[iC1].y-clusters2[iC2].y)/(clusters1[iC1].x-clusters2[iC2].x);
255	float r1 = Sqrt(clusters1[iC1].x * clusters1[iC1].x + clusters1[iC1].y * clusters1[iC1].y);
256	float r2 = Sqrt(clusters2[iC2].x * clusters2[iC2].x + clusters2[iC2].y * clusters2[iC2].y);
257	fDoublets[iL][fNDoublets[iL]].tanLambda = (clusters1[iC1].z-clusters2[iC2].z)/(r1-r2);
258	++fNDoublets[iL];
259	flag = false;
260	} else break;
261
262	}
263	if (flag) break;
264	}
265
266
267	// "Central" points
268	for ( int iC1 = 0 ; iC1 < fNClusters[iL] ; ++iC1 ) {
269
270	for ( int iC2 = 0; iC2 < fNClusters[iL+1] ; ++iC2 ) {
271
272	if( Abs( clusters1[iC1].phi - clusters2[iC2].phi ) < fPhiCut ) {
273	fDoublets[iL].push_back(nPlets(iC1,iC2));
274	fDoublets[iL][fNDoublets[iL]].tanPhi = (clusters1[iC1].y-clusters2[iC2].y)/(clusters1[iC1].x-clusters2[iC2].x);
275	float r1 = Sqrt(clusters1[iC1].x * clusters1[iC1].x + clusters1[iC1].y * clusters1[iC1].y);
276	float r2 = Sqrt(clusters2[iC2].x * clusters2[iC2].x + clusters2[iC2].y * clusters2[iC2].y);
277	fDoublets[iL][fNDoublets[iL]].tanLambda = (clusters1[iC1].z-clusters2[iC2].z)/(r1-r2);
278	++fNDoublets[iL];
279	} else if( clusters2[iC2].phi - clusters1[iC1].phi > fPhiCut ) break;
280
281	}
282
283	}
284
285	// 0 - 2Pi junction treatment (part II)
286	for ( int iC1 = fNClusters[iL]-1; iC1 > -1 ; --iC1 ) {
287	bool flag = true;
288
289	for ( int iC2 = 0; iC2 < fNClusters[iL+1] ; ++iC2 ) {
290
291	if( (TMath::TwoPi() - (clusters1[iC1].phi-clusters2[iC2].phi) ) < fPhiCut ) {
292	fDoublets[iL].push_back(nPlets(iC1,iC2));
293	fDoublets[iL][fNDoublets[iL]].tanPhi = (clusters1[iC1].y-clusters2[iC2].y)/(clusters1[iC1].x-clusters2[iC2].x);
294	float r1 = Sqrt(clusters1[iC1].x * clusters1[iC1].x + clusters1[iC1].y * clusters1[iC1].y);
295	float r2 = Sqrt(clusters2[iC2].x * clusters2[iC2].x + clusters2[iC2].y * clusters2[iC2].y);
296	fDoublets[iL][fNDoublets[iL]].tanLambda = (clusters1[iC1].z-clusters2[iC2].z)/(r1-r2);
297	++fNDoublets[iL];
298	flag = false;
299	} else break;
300
301	}
302
303	if (flag) break;
304	}
305
306	}
307
308	}