865c0bead16610964758ff7ac7e36c376b461427
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUVertexer.cxx
1 #include <Riostream.h>
2 #include <TBranch.h>
3 #include <TClonesArray.h>
4 #include <TMath.h>
5 #include <TString.h>
6 #include <TTree.h>
7 #include "AliESDVertex.h"
8 #include "AliITSUClusterLines.h"
9 #include "AliITSUClusterPix.h"
10 #include "AliITSUVertexer.h"
11 #include "AliLog.h"
12 #include "AliStrLine.h"
13 #include "AliVertexerTracks.h"
14
15 using TMath::Abs;
16 using TMath::Sqrt;
17 using TMath::ATan2;
18 using TMath::TwoPi;
19 using TMath::BubbleLow;
20
21 using std::cout;
22 using std::endl;
23
24 //////////////////////////////////////////////////////////////////////
25 // This class is used to compute the position of all the primary    // 
26 // vertices in a single event using the upgraded ITS.               //
27 // Optimizations ongoing.                                           //
28 // Origin puccio@to.infn.it  Feb. 20 2014                           //
29 //////////////////////////////////////////////////////////////////////
30
31 //_____________________________________________________________________________________________
32 AliITSUVertexer::AliITSUVertexer(Double_t phicut, Double_t zcut, Double_t paircut,Double_t clcut, Int_t cclcut) : AliVertexer(),
33                                                                                                                   fClusterContribCut(cclcut),
34                                                                                                                   fClusterCut(clcut),
35                                                                                                                   fClusterIndex(),
36                                                                                                                   fClusterPhi(),
37                                                                                                                   fClusters(),
38                                                                                                                   fLines("AliStrLine",1000),
39                                                                                                                   fLinesClusters("AliITSUClusterLines.h",1000),
40                                                                                                                   fLinesPhi(0),
41                                                                                                                   fNoClusters(0),
42                                                                                                                   fNoLines(0),
43                                                                                                                   fNoVertices(0),
44                                                                                                                   fPairCut(paircut),
45                                                                                                                   fPhiCut(phicut),
46                                                                                                                   fZCut(zcut),
47                                                                                                                   fUsedClusters(),
48                                                                                                                   fUsedLines(),
49                                                                                                                   fVertices(NULL) 
50 #ifdef MC_CHECK 
51 ,fGoodLines(0),fGoodLinesPhi(0),fParticleId(0)
52 #endif
53 {
54   // Standard I/O constructor
55 }
56
57 //_____________________________________________________________________________________________
58 AliITSUVertexer::~AliITSUVertexer() {
59   // Destructor
60   Reset();
61 }
62
63 //_____________________________________________________________________________________________
64 void AliITSUVertexer::FindVerticesForCurrentEvent() {
65   // Try to find all the primary vertices in the current 
66   fNoVertices=0;
67   FindTracklets();
68   if(fNoLines<2) { 
69     //fVertices.push_back(AliESDVertex());
70     return;// AliESDVertex();
71   }
72   
73   //  fVertices.push_back(AliVertexerTracks::TrackletVertexFinder(&fLines,1));
74   //fNoVertices=1;
75   fUsedLines=new Short_t[fNoLines];
76   for(UInt_t i=0;i<fNoLines;++i) fUsedLines[i]=-1;
77   
78   fNoClusters=0;
79   for(UInt_t i1=0;i1<fNoLines;++i1) {
80     if(fUsedLines[i1]!=-1) continue;
81     AliStrLine* line1 = (AliStrLine*)fLines.At(i1);
82     for(UInt_t i2=i1+1;i2<fNoLines;++i2) {
83       if(fUsedLines[i2]!=-1) continue;
84       AliStrLine* line2 = (AliStrLine*)fLines.At(i2);
85       if(line1->GetDCA(line2)<=fPairCut) {
86         //cout << fNoClusters <<" " << i1 << " " << i2 << " ";
87         new(fLinesClusters[fNoClusters])AliITSUClusterLines(i1,line1,i2,line2);
88         AliITSUClusterLines* current=(AliITSUClusterLines*)fLinesClusters.At(fNoClusters);
89         Double_t p[3];
90         current->GetVertex(p);
91         if((p[0]*p[0]+p[1]*p[1])>=4) { // Beam pipe check
92           fLinesClusters.RemoveAt(fNoClusters);
93           fLinesClusters.Compress();
94           break;
95         }
96         fUsedLines[i1]=fNoClusters;
97         fUsedLines[i2]=fNoClusters;
98         for(UInt_t i3=0;i3<fNoLines;++i3) {
99           if(fUsedLines[i3]!=-1) continue;
100           AliStrLine *line3 = (AliStrLine*)fLines.At(i3);
101           //cout << p[0] << " " << p[1] << " " << p[2] << endl;
102           //line3->PrintStatus();
103           if(line3->GetDistFromPoint(p)<=fPairCut) {
104             //cout << i3 << " ";
105             current->Add(i3,line3);
106             fUsedLines[i3]=fNoClusters;
107             current->GetVertex(p);
108           }
109         }
110         ++fNoClusters;
111         //cout << endl;
112         break;
113       }
114     }
115   }
116   
117   fLinesClusters.Sort();
118
119   for(UInt_t i0=0;i0<fNoClusters; ++i0) {
120     Double_t p0[3],p1[3];
121     AliITSUClusterLines *clu0 = (AliITSUClusterLines*)fLinesClusters.At(i0);
122     clu0->GetVertex(p0);
123     for(UInt_t i1=i0+1;i1<fNoClusters; ++i1) {
124       AliITSUClusterLines *clu1 = (AliITSUClusterLines*)fLinesClusters.At(i1);
125       clu1->GetVertex(p1);
126       if (TMath::Abs(p0[2]-p1[2])<=fClusterCut) {
127         Double_t distance=(p0[0]-p1[0])*(p0[0]-p1[0])+(p0[1]-p1[1])*(p0[1]-p1[1])+(p0[2]-p1[2])*(p0[2]-p1[2]);
128         //Bool_t flag=kFALSE;
129         if(distance<=fPairCut*fPairCut) {
130           UInt_t n=0;
131           Int_t *labels=clu1->GetLabels(n);
132           for(UInt_t icl=0; icl<n; ++icl) clu0->Add(labels[icl],(AliStrLine*)fLines.At(labels[icl]));
133           clu0->GetVertex(p0);
134           //flag=kTRUE;
135         }
136         fLinesClusters.RemoveAt(i1);
137         fLinesClusters.Compress();
138         fNoClusters--;
139         i1--;
140         //if(flag) i1=10;
141       }
142     }
143   }
144
145   fVertices=new AliESDVertex[fNoClusters];
146   for(UInt_t i0=0; i0<fNoClusters; ++i0) {
147     AliITSUClusterLines *clu0 = (AliITSUClusterLines*)fLinesClusters.At(i0);
148     Int_t size=clu0->GetSize();
149     if(size<fClusterContribCut&&fNoClusters>1) {
150       fLinesClusters.RemoveAt(i0);
151       fLinesClusters.Compress();
152       fNoClusters--;
153       continue;
154     } 
155     Double_t p0[3],cov[6];
156     clu0->GetVertex(p0);
157     clu0->GetCovMatrix(cov);
158     if((p0[0]*p0[0]+p0[1]*p0[1])<1.98*1.98) {
159       fVertices[fNoVertices++]=AliESDVertex(p0,cov,99999.,size);   
160     }
161   }
162   
163   return;// AliVertexerTracks::TrackletVertexFinder(&fLines,0);
164 }
165
166 //______________________________________________________________________
167 AliESDVertex* AliITSUVertexer::FindVertexForCurrentEvent(TTree *cluTree)
168 {
169 // Reconstruction of all the primary vertices in the event. It returns the vertex with the highest number of contributors.  
170   Reset();
171   for(Int_t i=0;i<3;++i) {
172     TBranch* br = cluTree->GetBranch(Form("ITSRecPoints%d",i));
173     if (!br) return NULL;
174     br->SetAddress(&fClusters[i]);
175   }    
176   cluTree->GetEntry(0);
177   
178   SortClusters();
179
180   FindVerticesForCurrentEvent();
181   if(fNoVertices<1) return NULL;
182   return new AliESDVertex(fVertices[0]);
183 }  
184
185 //_____________________________________________________________________________________________
186 Int_t AliITSUVertexer::MatchPoints(UShort_t layer, Double_t anchor, Double_t *p0, Double_t *p1) {
187   // Method for matching clusters with similar phi and z inside the range [zmin,zmax]
188   AliDebug(3,Form("Matching points on layer %i...\n",layer));
189
190   Double_t xyz[3][3];        // {{0: 'min', 1: 'mid', 2: 'max'},{0: 'x', 1: 'y', 2: 'z'}}
191   AliITSUClusterPix* cl[3];  // {0: 'min', 1: 'mid', 2: 'max'}
192   Double_t phi[3];           // {0: 'min', 1: 'mid', 2: 'max'}
193
194   Int_t nocl=fClusters[layer]->GetEntriesFast();
195   Int_t imin=0,imax=nocl!=0 ? nocl-1 : nocl,imid=(imax+imin)/2;
196
197   Double_t a=0,b=0,c=0;
198   if(layer==2) {
199     a=(p0[0]-p1[0])*(p0[0]-p1[0])+(p0[1]-p1[1])*(p0[1]-p1[1]);
200     b=(p1[0]-p0[0])*p0[0]+(p1[1]-p0[1])*p0[1];
201     c=p0[0]*p0[0]+p0[1]*p0[1];
202   }
203
204   Int_t flag=-1;
205
206   while (imax > imin) {
207     while(fUsedClusters[layer][fClusterIndex[layer][imin]]==kTRUE&&imin<imax) ++imin;
208     while(fUsedClusters[layer][fClusterIndex[layer][imax]]==kTRUE&&imax>0) --imax;
209     if(imax<imin) return -1;    
210     
211     cl[0] = (AliITSUClusterPix*)fClusters[layer]->At(fClusterIndex[layer][imin]);
212     cl[0]->GetGlobalXYZ(xyz[0]);
213     phi[0] = fClusterPhi[layer][fClusterIndex[layer][imin]];
214     AliDebug(4,Form("Cluster min: %i, %i, %f, %i\n",imin,cl[0]->GetLabel(0),phi[0]-anchor,fUsedClusters[layer][fClusterIndex[layer][imin]]));
215     if((Abs(phi[0]-anchor)<=fPhiCut||Abs(Abs(phi[0]-anchor)-TMath::TwoPi())<=fPhiCut)) {
216       AliDebug(4,Form("Ok min: %i \n",imin));
217       if(layer==2) {
218         if(fUsedClusters[layer][fClusterIndex[layer][imin]]) {
219           flag=1;
220           break;
221         }
222         c-=(xyz[0][0]*xyz[0][0]+xyz[0][1]*xyz[0][1]);
223         Double_t z=p0[2]+(p1[2]-p0[2])*(-b+Sqrt(b*b-a*c))/a;
224         Double_t zmin=z-fZCut,zmax=z+fZCut;
225         if(zmin<=xyz[0][2]&&xyz[0][2]<=zmax) {
226           AliDebug(4,Form("Ok Z: %i \n",imin));
227           return fClusterIndex[layer][imin];
228         }
229         AliDebug(4,Form("No Z match: %i \n",imin));
230         c=p0[0]*p0[0]+p0[1]*p0[1];
231         flag=1;
232         break;
233       } else if(!fUsedClusters[layer][fClusterIndex[layer][imin]]) return fClusterIndex[layer][imin];
234     }
235     
236     cl[2] = (AliITSUClusterPix*)fClusters[layer]->At(fClusterIndex[layer][imax]);
237     cl[2]->GetGlobalXYZ(xyz[2]);
238     phi[2] = fClusterPhi[layer][fClusterIndex[layer][imax]];
239     AliDebug(4,Form("Cluster max: %i, %i, %f, %i\n",imax,cl[2]->GetLabel(0),phi[2]-anchor,fUsedClusters[layer][fClusterIndex[layer][imax]]));
240     if((Abs(phi[2]-anchor)<=fPhiCut||Abs(Abs(phi[2]-anchor)-TMath::TwoPi())<=fPhiCut)) {
241       AliDebug(4,Form("Ok max: %i \n",imax));
242       if(layer==2) {
243         if(fUsedClusters[layer][fClusterIndex[layer][imax]]) {
244           flag=3;
245           break;
246         }
247         c-=(xyz[2][0]*xyz[2][0]+xyz[2][1]*xyz[2][1]);
248         Double_t z=p0[2]+(p1[2]-p0[2])*(-b+Sqrt(b*b-a*c))/a;
249         Double_t zmin=z-fZCut,zmax=z+fZCut;
250         if(zmin<=xyz[2][2]&&xyz[2][2]<=zmax) {
251           AliDebug(4,Form("Ok Z: %i \n",imax));
252           return fClusterIndex[layer][imax];
253         }
254         c=p0[0]*p0[0]+p0[1]*p0[1];
255         AliDebug(4,Form("No Z match: %i \n",imax));
256         flag=3;
257         break;
258       } else if(!fUsedClusters[layer][fClusterIndex[layer][imax]]) return fClusterIndex[layer][imax];
259     }
260     
261     imid=(imax+imin)/2;
262     if(imid==imin) return -1;
263     Int_t step=1,sign=-1;
264     while(fUsedClusters[layer][fClusterIndex[layer][imid]]) {
265       imid=imid+step;
266       sign*=-1;
267       step=(step+sign)*(-1);
268       if(imid>=imax||imid<=imin) return -1;
269     }
270     
271     cl[1] = (AliITSUClusterPix*)fClusters[layer]->At(fClusterIndex[layer][imid]);
272     cl[1]->GetGlobalXYZ(xyz[1]);
273     phi[1] = fClusterPhi[layer][fClusterIndex[layer][imid]];
274     AliDebug(4,Form("Cluster mid: %i, %i, %f, %i \n",imid,cl[1]->GetLabel(0),phi[1]-anchor,fUsedClusters[layer][fClusterIndex[layer][imid]]));
275     //cout << imin << " " << imid << " " << imax << endl;
276     //cout << fClusterIndex[layer][imid] << endl;
277     if((Abs(phi[1]-anchor)<=fPhiCut||Abs(Abs(phi[1]-anchor)-TMath::TwoPi())<=fPhiCut)) {
278       AliDebug(4,Form("Ok mid: %i \n",imid));
279       if(layer==2) {
280         c-=(xyz[1][0]*xyz[1][0]+xyz[1][1]*xyz[1][1]);
281         Double_t z=p0[2]+(p1[2]-p0[2])*(-b+Sqrt(b*b-a*c))/a;
282         Double_t zmin=z-fZCut,zmax=z+fZCut;
283         if(zmin<=xyz[1][2]&&xyz[1][2]<=zmax) {
284           AliDebug(4,Form("Ok Z: %i \n",imid));
285           return fClusterIndex[layer][imid];
286         }
287         AliDebug(4,Form("No Z match: %i \n",imid));
288         c=p0[0]*p0[0]+p0[1]*p0[1];
289         flag=2;
290         break;
291       } else if(!fUsedClusters[layer][fClusterIndex[layer][imid]]) return fClusterIndex[layer][imid];
292     }
293     
294     // determine which subarray to search
295     if (phi[1] < anchor) {
296       // change min index to search upper subarray
297       AliDebug(4,"Case minor\n");
298       imin = imid + 1;
299       --imax;
300     } else if (phi[1] > anchor) {
301       AliDebug(4,"Case greater\n");
302       // change max index to search lower subarray
303       ++imin;
304       imax = imid - 1;
305     } else if(!fUsedClusters[layer][fClusterIndex[layer][imid]]) {
306       return fClusterIndex[layer][imid];
307     } else return -1;
308   }
309
310   if(flag>-1) {
311     AliDebug(4,"Flag issued, starting forward backward check\n");
312     Int_t start=imid;
313     switch(flag) {
314     case 1:
315       start=imin;
316       break;
317     case 2:
318       start=imid;
319       break;
320     case 3:
321       start=imax;
322       break;
323     }
324    
325     Int_t curr=start-1;
326     Bool_t lap=kFALSE;
327     while(1){
328       if(curr==-1&&!lap) {
329         lap=kTRUE;
330         curr=nocl-1;
331       } else if(lap) break;
332       cl[0] = (AliITSUClusterPix*)fClusters[layer]->At(fClusterIndex[layer][curr]);
333       cl[0]->GetGlobalXYZ(xyz[0]);
334       phi[0] = fClusterPhi[layer][fClusterIndex[layer][curr]];
335       AliDebug(4,Form("Looking backward: %i, %i, %f, %i\n",curr,cl[0]->GetLabel(0),phi[0]-anchor,fUsedClusters[layer][fClusterIndex[layer][curr]]));
336       if(!(Abs(phi[0]-anchor)<=fPhiCut||Abs(Abs(phi[0]-anchor)-TMath::TwoPi())<=fPhiCut)) break;
337       if(fUsedClusters[layer][fClusterIndex[layer][curr]]) {
338         curr--;
339         continue;
340       }
341       
342       AliDebug(4,Form("Ok backward: %i \n",curr));
343       
344       c-=(xyz[0][0]*xyz[0][0]+xyz[0][1]*xyz[0][1]);
345       Double_t z=p0[2]+(p1[2]-p0[2])*(-b+Sqrt(b*b-a*c))/a;
346       Double_t zmin=z-fZCut,zmax=z+fZCut;
347       if(zmin<=xyz[0][2]&&xyz[0][2]<=zmax) {
348         AliDebug(4,Form("Ok Z: %i \n",curr));
349         return fClusterIndex[layer][curr];
350       }
351       AliDebug(4,Form("No Z match: %i \n",curr));
352       c=p0[0]*p0[0]+p0[1]*p0[1];
353       curr--;
354     }
355     
356     lap=kFALSE;
357     curr=start+1;
358     while(1){
359       if(curr==nocl&&!lap) {
360         curr=0;
361         lap=kTRUE;
362       } else if(lap) break;
363       cl[0] = (AliITSUClusterPix*)fClusters[layer]->At(fClusterIndex[layer][curr]);
364       cl[0]->GetGlobalXYZ(xyz[0]);
365       phi[0] = fClusterPhi[layer][fClusterIndex[layer][curr]];
366       AliDebug(4,Form("Looking forward: %i, %i, %f, %i\n",curr,cl[0]->GetLabel(0),phi[0]-anchor,fUsedClusters[layer][fClusterIndex[layer][curr]]));
367       if(!(Abs(phi[0]-anchor)<=fPhiCut||Abs(Abs(phi[0]-anchor)-TMath::TwoPi())<=fPhiCut)) break;
368       if(fUsedClusters[layer][fClusterIndex[layer][curr]]) {
369         curr++;
370         continue;
371       }
372       AliDebug(4,Form("Ok forward: %i \n",curr));
373       c-=(xyz[0][0]*xyz[0][0]+xyz[0][1]*xyz[0][1]);
374       Double_t z=p0[2]+(p1[2]-p0[2])*(-b+Sqrt(b*b-a*c))/a;
375       Double_t zmin=z-fZCut,zmax=z+fZCut;
376       if(zmin<=xyz[0][2]&&xyz[0][2]<=zmax) {
377         AliDebug(4,Form("Ok Z: %i \n",curr));
378         return fClusterIndex[layer][curr];
379       }
380       AliDebug(4,Form("No Z match: %i \n",curr));
381       c=p0[0]*p0[0]+p0[1]*p0[1];
382       curr++;
383     }
384   }
385   
386   if(imax==imin&&imax!=0) {
387     cl[0] = (AliITSUClusterPix*)fClusters[layer]->At(fClusterIndex[layer][imin]);
388     cl[0]->GetGlobalXYZ(xyz[0]);
389     phi[0] = fClusterPhi[layer][fClusterIndex[layer][imin]];
390     AliDebug(4,Form("Cluster eq: %i, %i, %f, %i\n",imin,cl[0]->GetLabel(0),phi[0]-anchor,fUsedClusters[layer][fClusterIndex[layer][imin]]));
391     if((Abs(phi[0]-anchor)<=fPhiCut||Abs(Abs(phi[0]-anchor)-TMath::TwoPi())<=fPhiCut)&&!fUsedClusters[layer][fClusterIndex[layer][imin]]) {
392       AliDebug(4,Form("Ok eq: %i \n",imin));
393       if(layer==2) {
394         c-=(xyz[0][0]*xyz[0][0]+xyz[0][1]*xyz[0][1]);
395         Double_t z=p0[2]+(p1[2]-p0[2])*(-b+Sqrt(b*b-a*c))/a;
396         Double_t zmin=z-fZCut,zmax=z+fZCut;
397         if(zmin<=xyz[0][2]&&xyz[0][2]<=zmax) {
398           AliDebug(4,Form("Ok Z eq: %i \n",imin));
399           return fClusterIndex[layer][imin];
400         }       
401         AliDebug(4,Form("No Z eq: %i \n",imin));
402         c=p0[0]*p0[0]+p0[1]*p0[1];
403       } else return fClusterIndex[layer][imin];
404     }
405   }
406   // no match found :(
407   return -1;
408
409 }
410
411 //_____________________________________________________________________________________________
412 void AliITSUVertexer::PrintStatus() const {
413   // Prints all the cuts and important data members for the current status
414   cout << "Cut on phi: " << fPhiCut << endl;
415   cout << "Cut on z: " << fZCut << endl;
416 }
417
418 //_____________________________________________________________________________________________
419 void AliITSUVertexer::Reset() {
420   // Resets the vertexer for a new event (or for its destruction) 
421   AliDebug(2,"Resetting the vertexer...\n");
422   fNoVertices=0;
423   for(Int_t i=0;i<3;++i) {
424     delete[] fUsedClusters[i];
425     delete[] fClusterIndex[i];
426     delete[] fClusterPhi[i];
427   }
428   if(fNoLines>2) {
429     delete []fUsedLines;
430   }
431
432   delete[] fVertices;
433
434   fLinesPhi=0;
435   fLines.Clear();
436   fLinesClusters.Clear();
437   
438   #ifdef MC_CHECK
439   fGoodLines=0;
440   fGoodLinesPhi=0;
441   delete[] fParticleId;
442   #endif
443
444   //delete fVertices;
445 }
446
447 //_____________________________________________________________________________________________
448 void AliITSUVertexer::FindTracklets() {
449   // It combines recpoints over the first three layers to define a list of tracklets
450   // Here it uses the ordered list of points. Second and third layers are ordered using phi
451   AliDebug(2,"Calling the trackleter...\n");
452
453   UInt_t noPntL[3]; 
454   for(UShort_t i=0;i<3;++i) {
455     noPntL[i]=fClusters[i]->GetEntries();
456     fUsedClusters[i]=new Bool_t[noPntL[i]];
457     for(UInt_t ii=0;ii<noPntL[i];++ii) fUsedClusters[i][ii]=kFALSE;
458   }
459   
460   #ifdef MC_CHECK
461   fParticleId=new UInt_t[noPntL[0]];
462   #endif
463
464   fNoLines=0;
465   UInt_t nolinesphi=0;
466   Double_t p0[3],p1[3],pp2[3];
467   for(UInt_t i0=0;i0<noPntL[0];++i0) {
468     if(fUsedClusters[0][i0]) continue;
469     AliITSUClusterPix* cluster0=(AliITSUClusterPix*)fClusters[0]->At(i0);
470     cluster0->GetGlobalXYZ(p0);
471     vector<Int_t> tmp;
472     Int_t i1=0;
473     Int_t label0=cluster0->GetLabel(0);  
474     AliDebug(4,Form("Layer 0: %i\n",label0));
475     while(i1>=0) {
476       i1 = MatchPoints(1,fClusterPhi[0][i0]);
477       if(i1<0) break;
478       
479       ++nolinesphi;
480       AliITSUClusterPix* cluster1=(AliITSUClusterPix*)fClusters[1]->At(i1);
481       cluster1->GetGlobalXYZ(p1);
482       tmp.push_back(i1);
483       fUsedClusters[1][i1]=kTRUE;
484
485       Int_t i2 = MatchPoints(2,fClusterPhi[1][i1],p0,p1);
486       if(i2<0) continue;
487       
488       #ifdef MC_CHECK
489       CheckMC(i0,i1,i2);
490       #endif
491
492       AliITSUClusterPix* cluster2=(AliITSUClusterPix*)fClusters[2]->At(i2);
493       cluster2->GetGlobalXYZ(pp2);
494
495       fUsedClusters[0][i0]=kTRUE;
496       fUsedClusters[2][i2]=kTRUE;
497       
498       // Errors to be checked... 
499       
500       Float_t cov0[6],cov1[6],cov2[6];
501       cluster0->GetGlobalCov(cov0);
502       cluster1->GetGlobalCov(cov1);
503       cluster2->GetGlobalCov(cov2);
504       
505       //Error on tracklet direction near the vertex
506       Double_t rad1=TMath::Sqrt(p0[0]*p0[0]+p0[1]*p0[1]);
507       Double_t rad2=TMath::Sqrt(p1[0]*p1[0]+p1[1]*p1[1]);
508       Double_t factor=(rad1+rad2)/(rad2-rad1);
509       
510       //Curvature error
511       Double_t curvErr=0;
512       Double_t bField=0.5;
513       Double_t meanPtSelTrk=0.630;
514       Double_t curvRadius=meanPtSelTrk/(0.3*bField)*100; //cm 
515       Double_t dRad=TMath::Sqrt((p0[0]-p1[0])*(p0[0]-p1[0])+(p0[1]-p1[1])*(p0[1]-p1[1]));
516       Double_t aux=dRad/2.+rad1;
517       curvErr=TMath::Sqrt(curvRadius*curvRadius-dRad*dRad/4.)-TMath::Sqrt(curvRadius*curvRadius-aux*aux); //cm
518       
519       Double_t sq[3],wmat[9]={1,0,0,0,1,0,0,0,1};
520       sq[0]=(cov0[0]+curvErr*curvErr/2.)*factor*factor;//cov1[0]+cov2[0]);
521       sq[1]=(cov0[3]+curvErr*curvErr/2.)*factor*factor;//cov1[3]+cov2[3]);
522       sq[2]=(cov0[5])*factor*factor;//cov1[5]+cov2[5]);     
523
524       // Multiple scattering
525       Double_t meanPSelTrk=0.875;
526       Double_t pOverMass=meanPSelTrk/0.140;
527       Double_t beta2=pOverMass*pOverMass/(1+pOverMass*pOverMass);
528       Double_t p2=meanPSelTrk*meanPSelTrk;
529       Double_t rBP=1.98; // Beam pipe radius
530       Double_t dBP=0.08/35.3; // 800 um of Be
531       Double_t dL1=0.01; //approx. 1% of radiation length  
532       Double_t theta2BP=14.1*14.1/(beta2*p2*1e6)*dBP;
533       Double_t theta2L1=14.1*14.1/(beta2*p2*1e6)*dL1;
534       Double_t rtantheta1=(rad2-rad1)*TMath::Tan(TMath::Sqrt(theta2L1));
535       Double_t rtanthetaBP=(rad1-rBP)*TMath::Tan(TMath::Sqrt(theta2BP));
536       for(Int_t ico=0; ico<3;ico++){    
537         sq[ico]+=rtantheta1*rtantheta1*factor*factor/3.;
538         sq[ico]+=rtanthetaBP*rtanthetaBP*factor*factor/3.;
539       }
540       
541       if(sq[0]!=0) wmat[0]=1/sq[0];
542       if(sq[1]!=0) wmat[4]=1/sq[1];
543       if(sq[2]!=0) wmat[8]=1/sq[2];
544
545       /*Int_t label0=cluster0->GetLabel(0);
546       Int_t label1=cluster1->GetLabel(0);
547       Int_t label2=cluster2->GetLabel(0);*/
548       //cout << label0 << " " << label1 << " "<<label2 << endl;
549
550       //printf("%f\t%f\t%f\n-\t%f\t%f\n-\t-\t%f\n",wmat[0],wmat[1],wmat[2],wmat[3],wmat[4],wmat[5]);
551       //gSystem->Exec(Form("echo %i >> trkl",cluster0->GetLabel(0)+cluster1->GetLabel(0)));
552       
553       new(fLines[fNoLines++])AliStrLine(p0,sq,wmat,p1,kTRUE);
554       break;
555     }
556     for(UInt_t itmp=0; itmp<tmp.size(); ++itmp) fUsedClusters[1][tmp.at(itmp)]=kFALSE;
557     tmp.clear();
558     //((AliStrLine*)fLines[fNoLines-1])->PrintStatus();
559     //printf("(%f,%f,%f) and (%f,%f,%f)\n",p0[0],p0[1],p0[2],p1[0],p1[1],p1[2]);
560   }
561   
562   fLinesPhi=nolinesphi;
563   //cout << Form("=======================================================\nNolinesphi: %i, nolines: %i\nGood Nolinesphi: %i, good nolines: %i\n",nolinesphi,fNoLines,fGoodLinesPhi,fGoodLines);
564 }
565  
566 //___________________________________________________________________________
567 void AliITSUVertexer::SortClusters() { 
568   // Reading of the clusters on the first three layer of the upgraded ITS
569   for(Int_t i=0;i<3;++i) {
570     TClonesArray *clr=fClusters[i];
571     Int_t nocl=clr->GetEntriesFast();
572     if(nocl==0) {
573       fClusterPhi[i]=new Double_t[1];
574       fClusterPhi[i][0]=-999999;
575       fClusterIndex[i]=new Int_t[1];
576       fClusterIndex[i][0]=0;
577     } else {
578       fClusterPhi[i]=new Double_t[nocl];
579       fClusterIndex[i]=new Int_t[nocl];
580       for(Int_t k=0;k<nocl;++k) {
581         AliITSUClusterPix* cl=(AliITSUClusterPix*)clr->At(k);
582         Double_t pt[3]; 
583         cl->GetGlobalXYZ(pt);
584         fClusterPhi[i][k]=ATan2(pt[1],pt[0]);
585       }
586       BubbleLow(nocl,fClusterPhi[i],fClusterIndex[i]);
587     }
588   }
589 }
590
591 #ifdef MC_CHECK
592 //___________________________________________________________________________
593 Bool_t AliITSUVertexer::CheckMC(UInt_t i0, UInt_t i1, UInt_t i2) {
594   // Debugging function
595   //cout << "Checking MC truth" << endl;
596   Int_t label=0;
597   Bool_t flag=kFALSE;
598   AliITSUClusterPix* p0=(AliITSUClusterPix*)fClusters[0]->At(i0);
599   AliITSUClusterPix* p1=(AliITSUClusterPix*)fClusters[1]->At(i1);
600   for(Int_t i=0; i<3; ++i) {
601     label=p0->GetLabel(i);
602     for(Int_t j=0; j<3; ++j) {
603       if(label==p1->GetLabel(j)&&label>0) {
604         fGoodLinesPhi++;
605         flag=kTRUE;                           
606         break;
607       }
608     }
609     if(flag) break;
610   }
611  
612   if(!flag) return kFALSE;
613   
614   AliITSUClusterPix* p2=(AliITSUClusterPix*)fClusters[2]->At(i2);
615   for(Int_t j=0; j<3; ++j) {
616     if(label==p2->GetLabel(j)) {
617       fParticleId[fGoodLines]=label;
618       fGoodLines++;
619       //cout << label << endl;
620       return kTRUE;
621     }
622   }
623   return kFALSE;
624 }
625 #endif