]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/UPGRADE/AliITSUVertexer.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
minor fix
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUVertexer.cxx
CommitLineData
30da2335 1#include <Riostream.h>
aea3f062 2#include <TBranch.h>
3#include <TClonesArray.h>
30da2335 4#include <TMath.h>
5#include <TString.h>
aea3f062 6#include <TTree.h>
b824a5f5 7#include "AliESDVertex.h"
aea3f062 8#include "AliITSUClusterLines.h"
30da2335 9#include "AliITSUClusterPix.h"
10#include "AliITSUVertexer.h"
aea3f062 11#include "AliLog.h"
30da2335 12#include "AliStrLine.h"
13#include "AliVertexerTracks.h"
b824a5f5 14
30da2335 15using TMath::Abs;
16using TMath::Sqrt;
17using TMath::ATan2;
18using TMath::TwoPi;
19using TMath::BubbleLow;
b824a5f5 20
5e2539d9 21using std::cout;
22using std::endl;
23
30da2335 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//////////////////////////////////////////////////////////////////////
b824a5f5 30
30da2335 31//_____________________________________________________________________________________________
32AliITSUVertexer::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(),
b82baa0e 49 fVertices(NULL)
30da2335 50#ifdef MC_CHECK
51,fGoodLines(0),fGoodLinesPhi(0),fParticleId(0)
52#endif
b824a5f5 53{
30da2335 54 // Standard I/O constructor
b824a5f5 55}
56
30da2335 57//_____________________________________________________________________________________________
b824a5f5 58AliITSUVertexer::~AliITSUVertexer() {
59 // Destructor
30da2335 60 Reset();
61}
62
63//_____________________________________________________________________________________________
64void AliITSUVertexer::FindVerticesForCurrentEvent() {
65 // Try to find all the primary vertices in the current
66 fNoVertices=0;
67 FindTracklets();
68 if(fNoLines<2) {
b82baa0e 69 //fVertices.push_back(AliESDVertex());
30da2335 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 }
30da2335 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
b82baa0e 145 fVertices=new AliESDVertex[fNoClusters];
30da2335 146 for(Int_t i0=fNoClusters-1;i0>=0; --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) {
b82baa0e 159 fVertices[fNoVertices++]=AliESDVertex(p0,cov,99999.,size);
30da2335 160 }
161 }
30da2335 162
163 return;// AliVertexerTracks::TrackletVertexFinder(&fLines,0);
b824a5f5 164}
165
166//______________________________________________________________________
aea3f062 167AliESDVertex* AliITSUVertexer::FindVertexForCurrentEvent(TTree *cluTree)
b824a5f5 168{
aea3f062 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();
b82baa0e 181 if(fNoVertices<1) return NULL;
182 return new AliESDVertex(fVertices[0]);
b824a5f5 183}
184
30da2335 185//_____________________________________________________________________________________________
186Int_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//_____________________________________________________________________________________________
412void 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//_____________________________________________________________________________________________
419void 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
b82baa0e 432 delete[] fVertices;
433
30da2335 434 fLinesPhi=0;
435 fLines.Clear();
436 fLinesClusters.Clear();
30da2335 437
438 #ifdef MC_CHECK
439 fGoodLines=0;
440 fGoodLinesPhi=0;
441 delete[] fParticleId;
442 #endif
443
444 //delete fVertices;
445}
446
447//_____________________________________________________________________________________________
448void 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//___________________________________________________________________________
aea3f062 567void AliITSUVertexer::SortClusters() {
30da2335 568 // Reading of the clusters on the first three layer of the upgraded ITS
aea3f062 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]);
30da2335 587 }
30da2335 588 }
589}
590
591#ifdef MC_CHECK
592//___________________________________________________________________________
593Bool_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