]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - ITS/AliITSVertexer3D.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
load AddTaskEmcalJet.C
[u/mrichter/AliRoot.git] / ITS / AliITSVertexer3D.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 2006-2008, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15#include <TTree.h>
16#include "AliRunLoader.h"
17#include "AliESDVertex.h"
18#include "AliLog.h"
19#include "AliStrLine.h"
20#include "AliTracker.h"
21#include "AliITSDetTypeRec.h"
22#include "AliITSRecPoint.h"
23#include "AliITSRecPointContainer.h"
24#include "AliITSgeomTGeo.h"
25#include "AliVertexerTracks.h"
26#include "AliITSVertexer3D.h"
27#include "AliITSVertexerZ.h"
28#include "AliITSSortTrkl.h"
29/////////////////////////////////////////////////////////////////
30// this class implements a method to determine
31// the 3 coordinates of the primary vertex
32// optimized for
33// p-p collisions
34////////////////////////////////////////////////////////////////
35
36const Int_t AliITSVertexer3D::fgkMaxNumOfClDefault = 300;
37const Int_t AliITSVertexer3D::fgkMaxNumOfClRebinDefault = 500;
38const Int_t AliITSVertexer3D::fgkMaxNumOfClDownscaleDefault = 1000;
39const Float_t AliITSVertexer3D::fgk3DBinSizeDefault = 0.1;
40
41ClassImp(AliITSVertexer3D)
42
43/* $Id$ */
44
45//______________________________________________________________________
46AliITSVertexer3D::AliITSVertexer3D(Double_t zcut):
47 AliITSVertexer(),
48 fLines("AliStrLine",1000),
49 fVert3D(),
50 fCoarseDiffPhiCut(0.),
51 fFineDiffPhiCut(0.),
52 fCutOnPairs(0.),
53 fCoarseMaxRCut(0.),
54 fMaxRCut(0.),
55 fMaxRCut2(0.),
56 fZCutDiamond(0.),
57 fMaxZCut(0.),
58 fDCAcut(0.),
59 fDiffPhiMax(0.),
60 fMeanPSelTrk(0.),
61 fMeanPtSelTrk(0.),
62 fUsedCluster(kMaxCluPerMod*kNSPDMod),
63 fZHisto(0),
64 fDCAforPileup(0.),
65 fDiffPhiforPileup(0.),
66 fBinSizeR(0.),
67 fBinSizeZ(0.),
68 fPileupAlgo(0),
69 fMaxNumOfCl(fgkMaxNumOfClDefault),
70 fMaxNumOfClForRebin(fgkMaxNumOfClRebinDefault),
71 fMaxNumOfClForDownScale(fgkMaxNumOfClDownscaleDefault),
72 fNRecPLay1(0),
73 fNRecPLay2(0),
74 f3DBinSize(fgk3DBinSizeDefault),
75 fDoDownScale(kFALSE),
76 fGenerForDownScale(0),
77 f3DPeak(),
78 fHighMultAlgo(1),
79 fSwitchAlgorithm(kFALSE)
80{
81 // Default constructor
82 SetCoarseDiffPhiCut();
83 SetFineDiffPhiCut();
84 SetCutOnPairs();
85 SetCoarseMaxRCut();
86 SetMaxRCut();
87 SetMaxRCutAlgo2();
88 if(zcut>0.){
89 SetZCutDiamond(zcut);
90 }
91 else {
92 SetZCutDiamond();
93 }
94 SetMaxZCut();
95 SetDCACut();
96 SetDiffPhiMax();
97 SetMeanPSelTracks();
98 SetMeanPtSelTracks();
99 SetMinDCAforPileup();
100 SetDeltaPhiforPileup();
101 SetPileupAlgo();
102 SetBinSizeR();
103 SetBinSizeZ();
104 fGenerForDownScale=new TRandom3(987654321);
105}
106
107//______________________________________________________________________
108AliITSVertexer3D::AliITSVertexer3D(TRootIOCtor*):
109 AliITSVertexer(),
110 fLines("AliStrLine",1000),
111 fVert3D(),
112 fCoarseDiffPhiCut(0.),
113 fFineDiffPhiCut(0.),
114 fCutOnPairs(0.),
115 fCoarseMaxRCut(0.),
116 fMaxRCut(0.),
117 fMaxRCut2(0.),
118 fZCutDiamond(0.),
119 fMaxZCut(0.),
120 fDCAcut(0.),
121 fDiffPhiMax(0.),
122 fMeanPSelTrk(0.),
123 fMeanPtSelTrk(0.),
124 fUsedCluster(kMaxCluPerMod*kNSPDMod),
125 fZHisto(0),
126 fDCAforPileup(0.),
127 fDiffPhiforPileup(0.),
128 fBinSizeR(0.),
129 fBinSizeZ(0.),
130 fPileupAlgo(0),
131 fMaxNumOfCl(fgkMaxNumOfClDefault),
132 fMaxNumOfClForRebin(fgkMaxNumOfClRebinDefault),
133 fMaxNumOfClForDownScale(fgkMaxNumOfClDownscaleDefault),
134 fNRecPLay1(0),
135 fNRecPLay2(0),
136 f3DBinSize(fgk3DBinSizeDefault),
137 fDoDownScale(kFALSE),
138 fGenerForDownScale(0),
139 f3DPeak(),
140 fHighMultAlgo(1),
141 fSwitchAlgorithm(kFALSE)
142{
143 // I/O constructor
144
145
146}
147
148//______________________________________________________________________
149AliITSVertexer3D::~AliITSVertexer3D() {
150 // Destructor
151 fLines.Clear("C");
152 if(fZHisto) delete fZHisto;
153 if(fGenerForDownScale) delete fGenerForDownScale;
154}
155
156//______________________________________________________________________
157void AliITSVertexer3D::ResetVert3D(){
158 // Reset the fVert3D object and reset the used clusters
159 ResetVertex();
160 fVert3D.SetXv(0.);
161 fVert3D.SetYv(0.);
162 fVert3D.SetZv(0.);
163 fVert3D.SetDispersion(0.);
164 fVert3D.SetNContributors(0);
165 fUsedCluster.ResetAllBits(0);
166}
167//______________________________________________________________________
168AliESDVertex* AliITSVertexer3D::FindVertexForCurrentEvent(TTree *itsClusterTree){
169 // Defines the AliESDVertex for the current event
170 ResetVert3D();
171 AliDebug(1,"FindVertexForCurrentEvent - 3D - PROCESSING NEXT EVENT");
172 fLines.Clear("C");
173 fCurrentVertex = NULL;
174
175 Int_t nolines = FindTracklets(itsClusterTree,0);
176 Int_t rc;
177 if(nolines>=2){
178 if(fSwitchAlgorithm) {
179 rc = Prepare3DVertexPbPb();
180 FindVertex3D(itsClusterTree);
181 } else {
182 rc=Prepare3DVertex(0);
183 if(fVert3D.GetNContributors()>0){
184 fLines.Clear("C");
185 nolines = FindTracklets(itsClusterTree,1);
186 if(nolines>=2){
187 rc=Prepare3DVertex(1);
188 if(fPileupAlgo == 2 && rc == 0) FindVertex3DIterative();
189 else if(fPileupAlgo!=2 && rc == 0) FindVertex3D(itsClusterTree);
190 if(rc!=0) fVert3D.SetNContributors(0); // exclude this vertex
191 }
192 }
193 }
194 }
195
196 if(!fCurrentVertex){
197 AliITSVertexerZ vertz(GetNominalPos()[0],GetNominalPos()[1]);
198 vertz.SetDetTypeRec(GetDetTypeRec());
199 AliDebug(1,"Call Vertexer Z\n");
200 vertz.SetLowLimit(-fZCutDiamond);
201 vertz.SetHighLimit(fZCutDiamond);
202 AliESDVertex* vtxz = vertz.FindVertexForCurrentEvent(itsClusterTree);
203 if(vtxz){
204 Double_t position[3]={GetNominalPos()[0],GetNominalPos()[1],vtxz->GetZv()};
205 Double_t covmatrix[6];
206 vtxz->GetCovMatrix(covmatrix);
207 Double_t chi2=99999.;
208 Int_t nContr=vtxz->GetNContributors();
209 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr);
210 fCurrentVertex->SetDispersion(vtxz->GetDispersion());
211 fCurrentVertex->SetTitle("vertexer: Z");
212 fCurrentVertex->SetName("SPDVertexZ");
213 delete vtxz;
214 }
215
216 }
217 if(fComputeMultiplicity) FindMultiplicity(itsClusterTree);
218 return fCurrentVertex;
219}
220
221//______________________________________________________________________
222void AliITSVertexer3D::FindVertex3D(TTree *itsClusterTree){
223 // Instantiates the fCurrentVertex object. calle by FindVertexForCurrenEvent
224 Double_t vRadius=TMath::Sqrt(fVert3D.GetXv()*fVert3D.GetXv()+fVert3D.GetYv()*fVert3D.GetYv());
225 if(vRadius<GetPipeRadius() && fVert3D.GetNContributors()>0){
226 Double_t position[3]={fVert3D.GetXv(),fVert3D.GetYv(),fVert3D.GetZv()};
227 Double_t covmatrix[6];
228 fVert3D.GetCovMatrix(covmatrix);
229 Double_t chi2=99999.;
230 Int_t nContr=fVert3D.GetNContributors();
231 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr);
232 fCurrentVertex->SetTitle("vertexer: 3D");
233 fCurrentVertex->SetName("SPDVertex3D");
234 fCurrentVertex->SetDispersion(fVert3D.GetDispersion());
235 fNoVertices=1;
236
237 switch(fPileupAlgo){
238 case 0: PileupFromZ(); break;
239 case 1: FindOther3DVertices(itsClusterTree); break;
240 case 3: break; // no pileup algo
241 default: AliError("Wrong pileup algorithm"); break;
242 }
243 if(fNoVertices==1){
244 delete[] fVertArray;
245 fVertArray = new AliESDVertex[1];
246 fVertArray[0]=(*fCurrentVertex);
247 }
248 }
249}
250
251//______________________________________________________________________
252void AliITSVertexer3D::FindVertex3DIterative(){
253 // find vertex if fPileupAlgo == 2
254
255 Int_t nLines=fLines.GetEntriesFast();
256 Int_t maxPoints=nLines*(nLines-1)/2;
257 Double_t* xP=new Double_t[maxPoints];
258 Double_t* yP=new Double_t[maxPoints];
259 Double_t* zP=new Double_t[maxPoints];
260 Int_t* index1=new Int_t[maxPoints];
261 Int_t* index2=new Int_t[maxPoints];
262 Double_t xbeam=fVert3D.GetXv();
263 Double_t ybeam=fVert3D.GetYv();
264
265 Int_t iPoint=0;
266 for(Int_t ilin1=0; ilin1<nLines; ilin1++){
267 AliStrLine *l1 = (AliStrLine*)fLines.At(ilin1);
268 for(Int_t ilin2=ilin1+1; ilin2<nLines; ilin2++){
269 AliStrLine *l2 = (AliStrLine*)fLines.At(ilin2);
270 Double_t dca=l1->GetDCA(l2);
271 if(dca > fDCAcut || dca<0.00001) continue;
272 Double_t point[3];
273 Int_t retc = l1->Cross(l2,point);
274 if(retc<0)continue;
275 Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
276 if(rad>fCoarseMaxRCut)continue;
277 Double_t distFromBeam=TMath::Sqrt((point[0]-xbeam)*(point[0]-xbeam)+(point[1]-ybeam)*(point[1]-ybeam));
278 if(distFromBeam>fMaxRCut2) continue;
279 xP[iPoint]=point[0];
280 yP[iPoint]=point[1];
281 zP[iPoint]=point[2];
282 index1[iPoint]=ilin1;
283 index2[iPoint]=ilin2;
284 iPoint++;
285 }
286 }
287 Int_t npoints=iPoint++;
288 Int_t index=0;
289 Short_t* mask=new Short_t[npoints];
290 for(Int_t ip=0;ip<npoints;ip++) mask[ip]=-1;
291
292 for(Int_t ip1=0;ip1<npoints;ip1++){
293 if(mask[ip1]==-1) mask[ip1]=index++;
294 for(Int_t ip2=ip1+1; ip2<npoints; ip2++){
295 if(mask[ip2]==mask[ip1] && mask[ip2]!=-1) continue;
296 Double_t dist2=(xP[ip1]-xP[ip2])*(xP[ip1]-xP[ip2]);
297 dist2+=(yP[ip1]-yP[ip2])*(yP[ip1]-yP[ip2]);
298 dist2+=(zP[ip1]-zP[ip2])*(zP[ip1]-zP[ip2]);
299 if(dist2<fCutOnPairs*fCutOnPairs){
300 if(mask[ip2]==-1) mask[ip2]=mask[ip1];
301 else{
302 for(Int_t ip=0; ip<npoints;ip++){
303 if(mask[ip]==mask[ip2]) mask[ip]=mask[ip1];
304 }
305 }
306 }
307 }
308 }
309
310
311 // Count multiplicity of trackelts in clusters
312 UInt_t* isIndUsed=new UInt_t[index+1];
313 for(Int_t ind=0;ind<index+1;ind++) isIndUsed[ind]=0;
314 for(Int_t ip=0; ip<npoints;ip++){
315 Int_t ind=mask[ip];
316 isIndUsed[ind]++;
317 }
318
319 // Count clusters/vertices and sort according to multiplicity
320 Int_t nClusters=0;
321 Int_t* sortedIndex=new Int_t[index+1];
322 for(Int_t ind1=0;ind1<index+1;ind1++){
323 if(isIndUsed[ind1]<=1) isIndUsed[ind1]=0;
324 else nClusters++;
325 UInt_t cap=9999999;
326 if(ind1>0) cap=isIndUsed[sortedIndex[ind1-1]];
327 UInt_t bigger=0;
328 Int_t biggerindex=-1;
329 for(Int_t ind2=0;ind2<index+1;ind2++){
330 Bool_t use=kTRUE;
331 for(Int_t ind3=0; ind3<ind1; ind3++)
332 if(ind2==sortedIndex[ind3]) use=kFALSE;
333 if(use && isIndUsed[ind2]>bigger && isIndUsed[ind2]<=cap){
334 bigger=isIndUsed[ind2];
335 biggerindex=ind2;
336 }
337 }
338 sortedIndex[ind1]=biggerindex;
339 }
340 AliDebug(3,Form("Number of clusters before merging = %d\n",nClusters));
341
342 // Assign lines to clusters/vertices and merge clusters which share 1 line
343 Int_t nClustersAfterMerge=nClusters;
344 Int_t* belongsTo=new Int_t[nLines];
345 for(Int_t ilin=0; ilin<nLines; ilin++) belongsTo[ilin]=-1;
346 for(Int_t iclu=0;iclu<nClusters;iclu++){
347 Int_t actualCluIndex=iclu;
348 for(Int_t ip=0; ip<npoints;ip++){
349 if(mask[ip]==sortedIndex[iclu]){
350 Int_t ind1=index1[ip];
351 if(belongsTo[ind1]==-1) belongsTo[ind1]=actualCluIndex;
352 else if(belongsTo[ind1]<actualCluIndex){
353 Int_t newCluIndex=belongsTo[ind1];
354 for(Int_t ilin=0; ilin<nLines; ilin++){
355 if(belongsTo[ilin]==actualCluIndex) belongsTo[ilin]=newCluIndex;
356 }
357 AliDebug(10,Form("Merged cluster %d with %d\n",actualCluIndex,newCluIndex));
358 actualCluIndex=newCluIndex;
359 nClustersAfterMerge--;
360 }
361 Int_t ind2=index2[ip];
362 if(belongsTo[ind2]==-1) belongsTo[ind2]=actualCluIndex;
363 else if(belongsTo[ind2]<actualCluIndex){
364 Int_t newCluIndex=belongsTo[ind2];
365 for(Int_t ilin=0; ilin<nLines; ilin++){
366 if(belongsTo[ilin]==actualCluIndex) belongsTo[ilin]=newCluIndex;
367 }
368 AliDebug(10,Form("Merged cluster %d with %d\n",actualCluIndex,newCluIndex));
369 actualCluIndex=newCluIndex;
370 nClustersAfterMerge--;
371 }
372 }
373 }
374 }
375 AliDebug(3,Form("Number of clusters after merging = %d\n",nClustersAfterMerge));
376
377 // Count lines associated to each cluster/vertex
378 UInt_t *cluSize=new UInt_t[nClusters];
379 for(Int_t iclu=0;iclu<nClusters;iclu++){
380 cluSize[iclu]=0;
381 for(Int_t ilin=0; ilin<nLines; ilin++){
382 if(belongsTo[ilin]==iclu) cluSize[iclu]++;
383 }
384 }
385
386 // Count good vertices (>1 associated tracklet)
387 UInt_t nGoodVert=0;
388 for(Int_t iclu=0;iclu<nClusters;iclu++){
389 AliDebug(3,Form("Vertex %d Size=%d\n",iclu,cluSize[iclu]));
390 if(cluSize[iclu]>1) nGoodVert++;
391 }
392
393 AliDebug(1,Form("Number of good vertices = %d\n",nGoodVert));
394 // Calculate vertex coordinates for each cluster
395 if(nGoodVert>0){
396 fVertArray = new AliESDVertex[nGoodVert];
397 Int_t iVert=0;
398 for(Int_t iclu=0;iclu<nClusters;iclu++){
399 Int_t size=cluSize[iclu];
400 if(size>1){
401 AliStrLine **arrlin = new AliStrLine*[size];
402 Int_t nFilled=0;
403 for(Int_t ilin=0; ilin<nLines; ilin++){
404 if(belongsTo[ilin]==iclu){
405 arrlin[nFilled++] = dynamic_cast<AliStrLine*>(fLines[ilin]);
406 }
407 }
408 AliDebug(3,Form("Vertex %d N associated tracklets = %d out of %d\n",iVert,size,nFilled));
409
410 fVertArray[iVert]=AliVertexerTracks::TrackletVertexFinder(arrlin,nFilled);
411 Double_t peak[3];
412 fVertArray[iVert].GetXYZ(peak);
413 AliStrLine **arrlin2 = new AliStrLine*[size];
414 Int_t nFilled2=0;
415 for(Int_t i=0; i<nFilled;i++){
416 AliStrLine *l1 = arrlin[i];
417 if(l1->GetDistFromPoint(peak)< fDCAcut)
418 arrlin2[nFilled2++] = dynamic_cast<AliStrLine*>(l1);
419 }
420 if(nFilled2>1){
421 AliDebug(3,Form("Vertex %d recalculated with %d tracklets\n",iVert,nFilled2));
422 fVertArray[iVert]=AliVertexerTracks::TrackletVertexFinder(arrlin2,nFilled2);
423 }
424 delete [] arrlin;
425 delete [] arrlin2;
426 ++iVert;
427 }
428 }
429
430 if(nGoodVert > 1){
431 fIsPileup = kTRUE;
432 fNTrpuv = fVertArray[1].GetNContributors();
433 fZpuv = fVertArray[1].GetZv();
434 }
435
436 Double_t vRadius=TMath::Sqrt(fVertArray[0].GetXv()*fVertArray[0].GetXv()+fVertArray[0].GetYv()*fVertArray[0].GetYv());
437 if(vRadius<GetPipeRadius() && fVertArray[0].GetNContributors()>0){
438 Double_t position[3]={fVertArray[0].GetXv(),fVertArray[0].GetYv(),fVertArray[0].GetZv()};
439 Double_t covmatrix[6];
440 fVertArray[0].GetCovMatrix(covmatrix);
441 Double_t chi2=99999.;
442 Int_t nContr=fVertArray[0].GetNContributors();
443 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr);
444 fCurrentVertex->SetTitle("vertexer: 3D");
445 fCurrentVertex->SetName("SPDVertex3D");
446 fCurrentVertex->SetDispersion(fVertArray[0].GetDispersion());
447 }
448 }
449
450 delete [] index1;
451 delete [] index2;
452 delete [] mask;
453 delete [] isIndUsed;
454 delete [] sortedIndex;
455 delete [] belongsTo;
456 delete [] cluSize;
457 delete [] xP;
458 delete [] yP;
459 delete [] zP;
460}
461//______________________________________________________________________
462void AliITSVertexer3D::FindVertex3DIterativeMM(){
463 // Defines the AliESDVertex for the current event
464 Int_t numsor=fLines.GetEntriesFast()*(fLines.GetEntriesFast()-1)/2;
465 //cout<<"AliITSVertexer3D::FindVertexForCurentEvent: Number of tracklets selected for vertexing "<<fLines.GetEntriesFast()<<"; Number of pairs: "<<numsor<<endl;
466 AliITSSortTrkl srt(fLines,numsor,fCutOnPairs,fCoarseMaxRCut);
467 srt.FindClusters();
468 AliInfo(Form("Number of vertices: %d",srt.GetNumberOfClusters()));
469
470 fNoVertices = srt.GetNumberOfClusters();
471 //printf("fNoVertices = %d \n",fNoVertices);
472 if(fNoVertices>0){
473 fVertArray = new AliESDVertex[fNoVertices];
474 for(Int_t kk=0; kk<srt.GetNumberOfClusters(); kk++){
475 Int_t size = 0;
476 Int_t *labels = srt.GetTrackletsLab(kk,size);
477 /*
478 Int_t *pairs = srt.GetClusters(kk);
479 Int_t nopai = srt.GetSizeOfCluster(kk);
480 cout<<"***** Vertex number "<<kk<<". Pairs: \n";
481 for(Int_t jj=0;jj<nopai;jj++){
482 cout<<pairs[jj]<<" - ";
483 if(jj>0 & jj%8==0)cout<<endl;
484 }
485 cout<<endl;
486 cout<<"***** Vertex number "<<kk<<". Labels: \n";
487 */
488 AliStrLine **tclo = new AliStrLine* [size];
489 for(Int_t jj=0;jj<size;jj++){
490 // cout<<labels[jj]<<" - ";
491 // if(jj>0 & jj%8==0)cout<<endl;
492 tclo[jj] = dynamic_cast<AliStrLine*>(fLines[labels[jj]]);
493 }
494 // cout<<endl;
495 delete []labels;
496 fVertArray[kk]=AliVertexerTracks::TrackletVertexFinder(tclo,size);
497 delete [] tclo;
498 // fVertArray[kk].PrintStatus();
499 if(kk == 1){
500 // at least one second vertex is present
501 fIsPileup = kTRUE;
502 fNTrpuv = fVertArray[kk].GetNContributors();
503 fZpuv = fVertArray[kk].GetZv();
504 }
505 }
506 Double_t vRadius=TMath::Sqrt(fVertArray[0].GetXv()*fVertArray[0].GetXv()+fVertArray[0].GetYv()*fVertArray[0].GetYv());
507 if(vRadius<GetPipeRadius() && fVertArray[0].GetNContributors()>0){
508 Double_t position[3]={fVertArray[0].GetXv(),fVertArray[0].GetYv(),fVertArray[0].GetZv()};
509 Double_t covmatrix[6];
510 fVertArray[0].GetCovMatrix(covmatrix);
511 Double_t chi2=99999.;
512 Int_t nContr=fVertArray[0].GetNContributors();
513 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr);
514 fCurrentVertex->SetTitle("vertexer: 3D");
515 fCurrentVertex->SetName("SPDVertex3D");
516 fCurrentVertex->SetDispersion(fVertArray[0].GetDispersion());
517 }
518 }
519
520}
521
522//______________________________________________________________________
523Bool_t AliITSVertexer3D::DistBetweenVertices(AliESDVertex &a, AliESDVertex &b, Double_t test, Double_t &dist){
524 // method to compare the distance between vertices a and b with "test"
525 //it returns kTRUE is the distance is less or equal to test
526 dist = (a.GetX()-b.GetX()) * (a.GetX()-b.GetX());
527 dist += (a.GetY()-b.GetY()) * (a.GetY()-b.GetY());
528 dist += (a.GetZ()-b.GetZ()) * (a.GetZ()-b.GetZ());
529 dist = TMath::Sqrt(dist);
530 if(dist <= test)return kTRUE;
531 return kFALSE;
532}
533
534
535//______________________________________________________________________
536Int_t AliITSVertexer3D::FindTracklets(TTree *itsClusterTree, Int_t optCuts){
537 // All the possible combinations between recpoints on layer 1and 2 are
538 // considered. Straight lines (=tracklets)are formed.
539 // The tracklets are processed in Prepare3DVertex
540
541 TClonesArray *itsRec = 0;
542 if(optCuts==0) fZHisto->Reset();
543 // gc1 are local and global coordinates for layer 1
544 Float_t gc1f[3]={0.,0.,0.};
545 Double_t gc1[3]={0.,0.,0.};
546 // gc2 are local and global coordinates for layer 2
547 Float_t gc2f[3]={0.,0.,0.};
548 Double_t gc2[3]={0.,0.,0.};
549 AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
550 rpcont->FetchClusters(0,itsClusterTree);
551 if(!rpcont->IsSPDActive()){
552 AliWarning("No SPD rec points found, 3D vertex not calculated");
553 return -1;
554 }
555
556 // Set values for cuts
557 Double_t xbeam=GetNominalPos()[0];
558 Double_t ybeam=GetNominalPos()[1];
559 Double_t zvert=0.;
560 Double_t deltaPhi=fCoarseDiffPhiCut;
561 Double_t deltaR=fCoarseMaxRCut;
562 Double_t dZmax=fZCutDiamond;
563 if(optCuts==1){
564 xbeam=fVert3D.GetXv();
565 ybeam=fVert3D.GetYv();
566 zvert=fVert3D.GetZv();
567 deltaPhi = fDiffPhiMax;
568 deltaR=fMaxRCut;
569 dZmax=fMaxZCut;
570 if(fPileupAlgo == 2){
571 dZmax=fZCutDiamond;
572 deltaR=fMaxRCut2;
573 }
574 } else if(optCuts==2){
575 xbeam=fVert3D.GetXv();
576 ybeam=fVert3D.GetYv();
577 deltaPhi = fDiffPhiforPileup;
578 deltaR=fMaxRCut;
579 }
580
581 fNRecPLay1=rpcont->GetNClustersInLayerFast(1);
582 fNRecPLay2=rpcont->GetNClustersInLayerFast(2);
583 if(fNRecPLay1 == 0 || fNRecPLay2 == 0){
584 AliDebug(1,Form("No RecPoints in at least one SPD layer (%d %d)",fNRecPLay1,fNRecPLay2));
585 return -1;
586 }
587 AliDebug(1,Form("RecPoints on Layer 1,2 = %d, %d\n",fNRecPLay1,fNRecPLay2));
588 fDoDownScale=kFALSE;
589 fSwitchAlgorithm=kFALSE;
590
591 Float_t factDownScal=1.;
592 Int_t origLaddersOnLayer2=fLadOnLay2;
593
594 switch(fHighMultAlgo){
595 case 0:
596 if(fNRecPLay1>fMaxNumOfClForDownScale || fNRecPLay2>fMaxNumOfClForDownScale){
597 if(optCuts==2) return -1; // do not try to search for pileup
598 SetLaddersOnLayer2(2);
599 fDoDownScale=kTRUE;
600 factDownScal=(Float_t)fMaxNumOfClForDownScale*(Float_t)fMaxNumOfClForDownScale/(Float_t)fNRecPLay1/(Float_t)fNRecPLay2;
601 if(optCuts==1){
602 factDownScal*=(fCoarseDiffPhiCut/fDiffPhiMax)*10;
603 if(factDownScal>1.){
604 fDoDownScale=kFALSE;
605 SetLaddersOnLayer2(origLaddersOnLayer2);
606 }
607 }
608 if(fDoDownScale)AliDebug(1,Form("Too many recpoints on SPD(%d %d ), downscale by %f",fNRecPLay1,fNRecPLay2,factDownScal));
609 }
610 break;
611 case 1:
612 if(fNRecPLay1>fMaxNumOfCl || fNRecPLay2>fMaxNumOfCl) {
613 if(optCuts==2) return -1; // do not try to search for pileup
614 fSwitchAlgorithm=kTRUE;
615 }
616 break;
617 default: break; // no pileup algo
618 }
619
620 if(!fDoDownScale && !fSwitchAlgorithm){
621 if(fNRecPLay1>fMaxNumOfClForRebin || fNRecPLay2>fMaxNumOfClForRebin){
622 SetLaddersOnLayer2(2);
623 }
624 }
625
626 Double_t a[3]={xbeam,ybeam,0.};
627 Double_t b[3]={xbeam,ybeam,10.};
628 AliStrLine zeta(a,b,kTRUE);
629 static Double_t bField=TMath::Abs(AliTracker::GetBz()/10.); //T
630 SetMeanPPtSelTracks(bField);
631
632 Int_t nolines = 0;
633 // Loop on modules of layer 1
634 Int_t firstL1 = TMath::Max(0,AliITSgeomTGeo::GetModuleIndex(1,1,1));
635 Int_t lastL1 = AliITSgeomTGeo::GetModuleIndex(2,1,1)-1;
636 for(Int_t modul1= firstL1; modul1<=lastL1;modul1++){ // Loop on modules of layer 1
637 if(!fUseModule[modul1]) continue;
638
639 UShort_t ladder=modul1/4+1; // ladders are numbered starting from 1
640 TClonesArray *prpl1=rpcont->UncheckedGetClusters(modul1);
641 Int_t nrecp1 = prpl1->GetEntries();
642 for(Int_t j=0;j<nrecp1;j++){
643 if(j>kMaxCluPerMod) continue;
644 UShort_t idClu1=modul1*kMaxCluPerMod+j;
645 if(fUsedCluster.TestBitNumber(idClu1)) continue;
646 if(fDoDownScale && !fSwitchAlgorithm){
647 if(fGenerForDownScale->Rndm()>factDownScal) continue;
648 }
649 AliITSRecPoint *recp1 = (AliITSRecPoint*)prpl1->At(j);
650 recp1->GetGlobalXYZ(gc1f);
651 for(Int_t ico=0;ico<3;ico++)gc1[ico]=gc1f[ico];
652
653 Double_t phi1 = TMath::ATan2(gc1[1]-ybeam,gc1[0]-xbeam);
654 if(phi1<0)phi1=2*TMath::Pi()+phi1;
655 for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){
656 for(Int_t k=0;k<4;k++){
657 Int_t ladmod=fLadders[ladder-1]+ladl2;
658 if(ladmod>AliITSgeomTGeo::GetNLadders(2)) ladmod=ladmod-AliITSgeomTGeo::GetNLadders(2);
659 Int_t modul2=AliITSgeomTGeo::GetModuleIndex(2,ladmod,k+1);
660 if(modul2<0)continue;
661 if(!fUseModule[modul2]) continue;
662 itsRec=rpcont->UncheckedGetClusters(modul2);
663 Int_t nrecp2 = itsRec->GetEntries();
664 for(Int_t j2=0;j2<nrecp2;j2++){
665 if(j2>kMaxCluPerMod) continue;
666 UShort_t idClu2=modul2*kMaxCluPerMod+j2;
667 if(fUsedCluster.TestBitNumber(idClu2)) continue;
668
669 AliITSRecPoint *recp2 = (AliITSRecPoint*)itsRec->At(j2);
670 recp2->GetGlobalXYZ(gc2f);
671 for(Int_t ico=0;ico<3;ico++)gc2[ico]=gc2f[ico];
672 Double_t phi2 = TMath::ATan2(gc2[1]-ybeam,gc2[0]-xbeam);
673 if(phi2<0)phi2=2*TMath::Pi()+phi2;
674 Double_t diff = TMath::Abs(phi2-phi1);
675 if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
676 if(optCuts==0 && diff<fDiffPhiforPileup){
677 Double_t r1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]);
678 Double_t zc1=gc1[2];
679 Double_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
680 Double_t zc2=gc2[2];
681 Double_t zr0=(r2*zc1-r1*zc2)/(r2-r1); //Z @ null radius
682 fZHisto->Fill(zr0);
683 }
684 if(diff>deltaPhi)continue;
685 AliStrLine line(gc1,gc2,kTRUE);
686 Double_t cp[3];
687 Int_t retcode = line.Cross(&zeta,cp);
688 if(retcode<0)continue;
689 Double_t dca = line.GetDCA(&zeta);
690 if(dca<0.) continue;
691 if(dca>deltaR)continue;
692 Double_t deltaZ=cp[2]-zvert;
693 if(TMath::Abs(deltaZ)>dZmax)continue;
694
695
696 if(nolines == 0){
697 if(fLines.GetEntriesFast()>0)fLines.Clear("C");
698 }
699 Float_t cov[6];
700 recp2->GetGlobalCov(cov);
701
702
703 Double_t rad1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]);
704 Double_t rad2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
705 Double_t factor=(rad1+rad2)/(rad2-rad1); //factor to account for error on tracklet direction
706
707 Double_t curvErr=0;
708 if(bField>0.00001){
709 Double_t curvRadius=fMeanPtSelTrk/(0.3*bField)*100; //cm
710 Double_t dRad=TMath::Sqrt((gc1[0]-gc2[0])*(gc1[0]-gc2[0])+(gc1[1]-gc2[1])*(gc1[1]-gc2[1]));
711 Double_t aux=dRad/2.+rad1;
712 curvErr=TMath::Sqrt(curvRadius*curvRadius-dRad*dRad/4.)-TMath::Sqrt(curvRadius*curvRadius-aux*aux); //cm
713 }
714 Double_t sigmasq[3];
715 sigmasq[0]=(cov[0]+curvErr*curvErr/2.)*factor*factor;
716 sigmasq[1]=(cov[3]+curvErr*curvErr/2.)*factor*factor;
717 sigmasq[2]=cov[5]*factor*factor;
718
719 // Multiple scattering
720 Double_t pOverMass=fMeanPSelTrk/0.140;
721 Double_t beta2=pOverMass*pOverMass/(1+pOverMass*pOverMass);
722 Double_t p2=fMeanPSelTrk*fMeanPSelTrk;
723 Double_t rBP=GetPipeRadius();
724 Double_t dBP=0.08/35.3; // 800 um of Be
725 Double_t dL1=0.01; //approx. 1% of radiation length
726 Double_t theta2BP=14.1*14.1/(beta2*p2*1e6)*dBP;
727 Double_t theta2L1=14.1*14.1/(beta2*p2*1e6)*dL1;
728 Double_t rtantheta1=(rad2-rad1)*TMath::Tan(TMath::Sqrt(theta2L1));
729 Double_t rtanthetaBP=(rad1-rBP)*TMath::Tan(TMath::Sqrt(theta2BP));
730 for(Int_t ico=0; ico<3;ico++){
731 sigmasq[ico]+=rtantheta1*rtantheta1*factor*factor/3.;
732 sigmasq[ico]+=rtanthetaBP*rtanthetaBP*factor*factor/3.;
733 }
734 Double_t wmat[9]={1.,0.,0.,0.,1.,0.,0.,0.,1.};
735 if(sigmasq[0]!=0.) wmat[0]=1./sigmasq[0];
736 if(sigmasq[1]!=0.) wmat[4]=1./sigmasq[1];
737 if(sigmasq[2]!=0.) wmat[8]=1./sigmasq[2];
738 new(fLines[nolines++])AliStrLine(gc1,sigmasq,wmat,gc2,kTRUE,idClu1,idClu2);
739
740 }
741 }
742 }
743 }
744 }
745
746 SetLaddersOnLayer2(origLaddersOnLayer2);
747
748 if(nolines == 0)return -2;
749 return nolines;
750}
751
752//______________________________________________________________________
753Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){
754 // Finds the 3D vertex information using tracklets
755 Int_t retcode = -1;
756 Double_t xbeam=GetNominalPos()[0];
757 Double_t ybeam=GetNominalPos()[1];
758 Double_t zvert=0.;
759 Double_t deltaR=fCoarseMaxRCut;
760 Double_t dZmax=fZCutDiamond;
761 if(optCuts==1){
762 xbeam=fVert3D.GetXv();
763 ybeam=fVert3D.GetYv();
764 zvert=fVert3D.GetZv();
765 deltaR=fMaxRCut;
766 dZmax=fMaxZCut;
767 if(fPileupAlgo == 2){
768 dZmax=fZCutDiamond;
769 deltaR=fMaxRCut2;
770 }
771 }else if(optCuts==2){
772 xbeam=fVert3D.GetXv();
773 ybeam=fVert3D.GetYv();
774 deltaR=fMaxRCut;
775 }
776
777 Double_t origBinSizeR=fBinSizeR;
778 Double_t origBinSizeZ=fBinSizeZ;
779 Bool_t rebinned=kFALSE;
780 if(fDoDownScale){
781 SetBinSizeR(0.05);
782 SetBinSizeZ(0.05);
783 rebinned=kTRUE;
784 }else{
785 if(optCuts==0 && (fNRecPLay1>fMaxNumOfClForRebin || fNRecPLay2>fMaxNumOfClForRebin)){
786 SetBinSizeR(0.1);
787 SetBinSizeZ(0.2);
788 rebinned=kTRUE;
789 }
790 }
791 Double_t rl=-fCoarseMaxRCut;
792 Double_t rh=fCoarseMaxRCut;
793 Double_t zl=-fZCutDiamond;
794 Double_t zh=fZCutDiamond;
795 Int_t nbr=(Int_t)((rh-rl)/fBinSizeR+0.0001);
796 Int_t nbz=(Int_t)((zh-zl)/fBinSizeZ+0.0001);
797 Int_t nbrcs=(Int_t)((rh-rl)/(fBinSizeR*2.)+0.0001);
798 Int_t nbzcs=(Int_t)((zh-zl)/(fBinSizeZ*2.)+0.0001);
799
800 TH3F *h3d = new TH3F("h3d","xyz distribution",nbr,rl,rh,nbr,rl,rh,nbz,zl,zh);
801 TH3F *h3dcs = new TH3F("h3dcs","xyz distribution",nbrcs,rl,rh,nbrcs,rl,rh,nbzcs,zl,zh);
802
803 // cleanup of the TCLonesArray of tracklets (i.e. fakes are removed)
804 Int_t vsiz = fLines.GetEntriesFast();
805 Int_t *validate = new Int_t [vsiz];
806 for(Int_t i=0; i<vsiz;i++)validate[i]=0;
807 for(Int_t i=0; i<vsiz-1;i++){
808 AliStrLine *l1 = (AliStrLine*)fLines.At(i);
809 for(Int_t j=i+1;j<fLines.GetEntriesFast();j++){
810 AliStrLine *l2 = (AliStrLine*)fLines.At(j);
811 Double_t dca=l1->GetDCA(l2);
812 if(dca > fDCAcut || dca<0.00001) continue;
813 Double_t point[3];
814 Int_t retc = l1->Cross(l2,point);
815 if(retc<0)continue;
816 Double_t deltaZ=point[2]-zvert;
817 if(TMath::Abs(deltaZ)>dZmax)continue;
818 Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
819 if(rad>fCoarseMaxRCut)continue;
820 Double_t deltaX=point[0]-xbeam;
821 Double_t deltaY=point[1]-ybeam;
822 Double_t raddist=TMath::Sqrt(deltaX*deltaX+deltaY*deltaY);
823 if(raddist>deltaR)continue;
824 validate[i]=1;
825 validate[j]=1;
826 h3d->Fill(point[0],point[1],point[2]);
827 h3dcs->Fill(point[0],point[1],point[2]);
828 }
829 }
830
831 Int_t numbtracklets=0;
832 for(Int_t i=0; i<vsiz;i++)if(validate[i]>=1)numbtracklets++;
833 if(numbtracklets<2){
834 delete [] validate;
835 delete h3d;
836 delete h3dcs;
837 SetBinSizeR(origBinSizeR);
838 SetBinSizeZ(origBinSizeZ);
839 return retcode;
840 }
841
842 for(Int_t i=0; i<fLines.GetEntriesFast();i++){
843 if(validate[i]<1)fLines.RemoveAt(i);
844 }
845 fLines.Compress();
846 AliDebug(1,Form("Number of tracklets (after compress)%d ",fLines.GetEntriesFast()));
847 delete [] validate;
848
849 // Exit here if Pileup Algorithm 2 has been chosen during second loop
850 if(fPileupAlgo == 2 && optCuts==1){
851 delete h3d;
852 delete h3dcs;
853 SetBinSizeR(origBinSizeR);
854 SetBinSizeZ(origBinSizeZ);
855 return 0;
856 }
857
858 // Find peaks in histos
859
860 Double_t peak[3]={0.,0.,0.};
861 Int_t ntrkl,ntimes;
862 FindPeaks(h3d,peak,ntrkl,ntimes);
863 delete h3d;
864 Double_t binsizer=(rh-rl)/nbr;
865 Double_t binsizez=(zh-zl)/nbz;
866 if(optCuts==0 && (ntrkl<=2 || ntimes>1)){
867 ntrkl=0;
868 ntimes=0;
869 FindPeaks(h3dcs,peak,ntrkl,ntimes);
870 binsizer=(rh-rl)/nbrcs;
871 binsizez=(zh-zl)/nbzcs;
872 if(ntrkl==1 || ntimes>1){
873 delete h3dcs;
874 SetBinSizeR(origBinSizeR);
875 SetBinSizeZ(origBinSizeZ);
876 return retcode;
877 }
878 }
879 delete h3dcs;
880
881 Double_t bs=(binsizer+binsizez)/2.;
882 for(Int_t i=0; i<fLines.GetEntriesFast();i++){
883 AliStrLine *l1 = (AliStrLine*)fLines.At(i);
884 if(l1->GetDistFromPoint(peak)>2.5*bs)fLines.RemoveAt(i);
885 }
886 fLines.Compress();
887 AliDebug(1,Form("Number of tracklets (after 2nd compression) %d",fLines.GetEntriesFast()));
888
889 // Finer Histo in limited range in case of high mult.
890 if(rebinned){
891 SetBinSizeR(0.01);
892 SetBinSizeZ(0.01);
893 Double_t xl=peak[0]-0.3;
894 Double_t xh=peak[0]+0.3;
895 Double_t yl=peak[1]-0.3;
896 Double_t yh=peak[1]+0.3;
897 zl=peak[2]-0.5;
898 zh=peak[2]+0.5;
899 Int_t nbxfs=(Int_t)((xh-xl)/fBinSizeR+0.0001);
900 Int_t nbyfs=(Int_t)((yh-yl)/fBinSizeR+0.0001);
901 Int_t nbzfs=(Int_t)((zh-zl)/fBinSizeZ+0.0001);
902
903 TH3F *h3dfs = new TH3F("h3dfs","xyz distribution",nbxfs,xl,xh,nbyfs,yl,yh,nbzfs,zl,zh);
904 for(Int_t i=0; i<fLines.GetEntriesFast()-1;i++){
905 AliStrLine *l1 = (AliStrLine*)fLines.At(i);
906 for(Int_t j=i+1;j<fLines.GetEntriesFast();j++){
907 AliStrLine *l2 = (AliStrLine*)fLines.At(j);
908 Double_t dca=l1->GetDCA(l2);
909 if(dca > fDCAcut || dca<0.00001) continue;
910 Double_t point[3];
911 Int_t retc = l1->Cross(l2,point);
912 if(retc<0)continue;
913 Double_t deltaZ=point[2]-zvert;
914 if(TMath::Abs(deltaZ)>dZmax)continue;
915 Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
916 if(rad>fCoarseMaxRCut)continue;
917 Double_t deltaX=point[0]-xbeam;
918 Double_t deltaY=point[1]-ybeam;
919 Double_t raddist=TMath::Sqrt(deltaX*deltaX+deltaY*deltaY);
920 if(raddist>deltaR)continue;
921 h3dfs->Fill(point[0],point[1],point[2]);
922 }
923 }
924 ntrkl=0;
925 ntimes=0;
926
927 Double_t newpeak[3]={0.,0.,0.};
928 FindPeaks(h3dfs,newpeak,ntrkl,ntimes);
929 if(ntimes==1){
930 for(Int_t iCoo=0; iCoo<3; iCoo++) peak[iCoo]=newpeak[iCoo];
931 binsizer=fBinSizeR;
932 binsizez=fBinSizeZ;
933 }
934 delete h3dfs;
935 bs=(binsizer+binsizez)/2.;
936 for(Int_t i=0; i<fLines.GetEntriesFast();i++){
937 AliStrLine *l1 = (AliStrLine*)fLines.At(i);
938 if(l1->GetDistFromPoint(peak)>2.5*bs)fLines.RemoveAt(i);
939 }
940 fLines.Compress();
941 AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines.GetEntriesFast()));
942 }
943 SetBinSizeR(origBinSizeR);
944 SetBinSizeZ(origBinSizeZ);
945
946
947 // Second selection loop
948
949
950 if(fLines.GetEntriesFast()>1){
951 retcode=0;
952 // find a first candidate for the primary vertex
953 fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
954 // make a further selection on tracklets based on this first candidate
955 fVert3D.GetXYZ(peak);
956 AliDebug(1,Form("FIRST V candidate: %f ; %f ; %f",peak[0],peak[1],peak[2]));
957 Int_t *validate2 = new Int_t [fLines.GetEntriesFast()];
958 for(Int_t i=0; i<fLines.GetEntriesFast();i++) validate2[i]=1;
959 for(Int_t i=0; i<fLines.GetEntriesFast();i++){
960 if(validate2[i]==0) continue;
961 AliStrLine *l1 = (AliStrLine*)fLines.At(i);
962 if(l1->GetDistFromPoint(peak)> fDCAcut)fLines.RemoveAt(i);
963 if(optCuts==2){ // temporarily only for pileup
964 for(Int_t j=i+1; j<fLines.GetEntriesFast();j++){
965 AliStrLine *l2 = (AliStrLine*)fLines.At(j);
966 if(l1->GetDCA(l2)<0.00001){
967 Int_t delta1=(Int_t)l1->GetIdPoint(0)-(Int_t)l2->GetIdPoint(0);
968 Int_t delta2=(Int_t)l1->GetIdPoint(1)-(Int_t)l2->GetIdPoint(1);
969 Int_t deltamod1=(Int_t)l1->GetIdPoint(0)/kMaxCluPerMod
970 -(Int_t)l2->GetIdPoint(0)/kMaxCluPerMod;
971 Int_t deltamod2=(Int_t)l1->GetIdPoint(1)/kMaxCluPerMod
972 -(Int_t)l2->GetIdPoint(1)/kMaxCluPerMod;
973 // remove tracklets sharing a point
974 if( (delta1==0 && deltamod2==0) ||
975 (delta2==0 && deltamod1==0) ) validate2[j]=0;
976 }
977 }
978 }
979 }
980 for(Int_t i=0; i<fLines.GetEntriesFast();i++){
981 if(validate2[i]==0) fLines.RemoveAt(i);
982 }
983 delete [] validate2;
984 fLines.Compress();
985 AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines.GetEntriesFast()));
986 if(fLines.GetEntriesFast()>1){// this new tracklet selection is used
987 fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
988 }
989 }
990 return retcode;
991}
992
993//________________________________________________________
994Int_t AliITSVertexer3D::Prepare3DVertexPbPb(){
995 // Finds the 3D vertex information in Pb-Pb events using tracklets
996 AliDebug(1,"High multiplicity event.\n");
997
998 Int_t nxy=(Int_t)(2.*fCoarseMaxRCut/f3DBinSize);
999 Double_t xymi= -nxy*f3DBinSize/2.;
1000 Double_t xyma= nxy*f3DBinSize/2.;
1001 Int_t nz=(Int_t)(2.*fZCutDiamond/f3DBinSize);
1002 Double_t zmi=-nz*f3DBinSize/2.;
1003 Double_t zma=nz*f3DBinSize/2.;
1004 Int_t nolines=fLines.GetEntriesFast();
1005 TH3F *h3dv = new TH3F("h3dv","3d tracklets",nxy,xymi,xyma,nxy,xymi,xyma,nz,zmi,zma);
1006
1007 for(Int_t itra=0; itra<nolines; itra++){
1008 Double_t wei = GetFraction(itra);
1009 //printf("tracklet %d ) - weight %f \n",itra,wei);
1010 if(wei>1.e-6){
1011 AliStrLine *str=(AliStrLine*)fLines.At(itra);
1012 Double_t t1,t2;
1013 if(str->GetParamAtRadius(fCoarseMaxRCut,t1,t2)){
1014 do{
1015 Double_t punt[3];
1016 str->ComputePointAtT(t1,punt);
1017 h3dv->Fill(punt[0],punt[1],punt[2],wei);
1018 t1+=f3DBinSize/3.;
1019 } while(t1<t2);
1020 }
1021 }
1022 }
1023 Int_t noftrk,noftim;
1024 FindPeaks(h3dv,f3DPeak,noftrk,noftim); // arg: histo3d, peak, # of contrib., # of other peak with same magnitude
1025
1026
1027 // Remove all the tracklets which are not passing near peak
1028
1029 while(nolines--){
1030 AliStrLine *str=(AliStrLine*)fLines.At(nolines);
1031 Double_t dist = str->GetDistFromPoint(f3DPeak);
1032 if(dist>(2.*f3DBinSize)) fLines.RemoveAt(nolines);
1033 }
1034 fLines.Compress();
1035 nolines=fLines.GetEntriesFast();
1036
1037 delete h3dv;
1038
1039 Int_t *validate2 = new Int_t [fLines.GetEntriesFast()];
1040 for(Int_t i=0; i<fLines.GetEntriesFast();i++) validate2[i]=1;
1041 for(Int_t i=0; i<fLines.GetEntriesFast();i++){
1042 if(validate2[i]==0) continue;
1043 AliStrLine *l1 = (AliStrLine*)fLines.At(i);
1044 if(l1->GetDistFromPoint(f3DPeak)> fDCAcut)fLines.RemoveAt(i);
1045 for(Int_t j=i+1; j<fLines.GetEntriesFast();j++){
1046 AliStrLine *l2 = (AliStrLine*)fLines.At(j);
1047 if(l1->GetDCA(l2)<0.00001){
1048 Int_t delta1=(Int_t)l1->GetIdPoint(0)-(Int_t)l2->GetIdPoint(0);
1049 Int_t delta2=(Int_t)l1->GetIdPoint(1)-(Int_t)l2->GetIdPoint(1);
1050 Int_t deltamod1=(Int_t)l1->GetIdPoint(0)/kMaxCluPerMod
1051 -(Int_t)l2->GetIdPoint(0)/kMaxCluPerMod;
1052 Int_t deltamod2=(Int_t)l1->GetIdPoint(1)/kMaxCluPerMod
1053 -(Int_t)l2->GetIdPoint(1)/kMaxCluPerMod;
1054 // remove tracklets sharing a point
1055 if( (delta1==0 && deltamod2==0) ||
1056 (delta2==0 && deltamod1==0) ) validate2[j]=0;
1057
1058 }
1059 }
1060 }
1061 for(Int_t i=0; i<fLines.GetEntriesFast();i++){
1062 if(validate2[i]==0) fLines.RemoveAt(i);
1063 }
1064
1065 delete [] validate2;
1066 fLines.Compress();
1067
1068
1069 AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines.GetEntriesFast()));
1070
1071 fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
1072 fVert3D.GetXYZ(f3DPeak);
1073
1074 return 0;
1075}
1076
1077//________________________________________________________
1078void AliITSVertexer3D::SetMeanPPtSelTracks(Double_t fieldTesla){
1079 // Sets mean values of Pt based on the field
1080 // for P (used in multiple scattering) the most probable value is used
1081 if(TMath::Abs(fieldTesla-0.5)<0.01){
1082 SetMeanPSelTracks(0.375);
1083 SetMeanPtSelTracks(0.630);
1084 }else if(TMath::Abs(fieldTesla-0.4)<0.01){
1085 SetMeanPSelTracks(0.375);
1086 SetMeanPtSelTracks(0.580);
1087 }else if(TMath::Abs(fieldTesla-0.2)<0.01){
1088 SetMeanPSelTracks(0.375);
1089 SetMeanPtSelTracks(0.530);
1090 }else if(fieldTesla<0.00001){
1091 SetMeanPSelTracks(0.375);
1092 SetMeanPtSelTracks(0.230);
1093 }else{
1094 SetMeanPSelTracks();
1095 SetMeanPtSelTracks();
1096 }
1097}
1098
1099//________________________________________________________
1100void AliITSVertexer3D::SetZCutDiamond(Double_t zcut){
1101 // The fiducial region along Z is set. The TH1 object pointed by
1102 // fZHisto is created
1103 if(TMath::AreEqualAbs(zcut,fZCutDiamond,1.e-10))return;
1104 fZCutDiamond=zcut;
1105 if(fZHisto) delete fZHisto;
1106 Double_t binsize=0.02; // default 200 micron
1107 Int_t nbins=static_cast<Int_t>(1+2*fZCutDiamond/binsize);
1108 fZHisto=new TH1F("hz","",nbins,-fZCutDiamond,-fZCutDiamond+binsize*nbins);
1109 fZHisto->SetDirectory(0);
1110}
1111
1112//________________________________________________________
1113void AliITSVertexer3D::FindPeaks(TH3F* histo, Double_t *peak, Int_t &nOfTracklets, Int_t &nOfTimes){
1114 // Finds bin with max contents in 3D histo of tracket intersections
1115 TAxis *xax = histo->GetXaxis();
1116 TAxis *yax = histo->GetYaxis();
1117 TAxis *zax = histo->GetZaxis();
1118 peak[0]=0.;
1119 peak[1]=0.;
1120 peak[2]=0.;
1121 nOfTracklets = 0;
1122 nOfTimes=0;
1123 Int_t peakbin[3]={0,0,0};
1124 Int_t peak2bin[3]={-1,-1,-1};
1125 Int_t bc2=-1;
1126 for(Int_t i=xax->GetFirst();i<=xax->GetLast();i++){
1127 Double_t xval = xax->GetBinCenter(i);
1128 for(Int_t j=yax->GetFirst();j<=yax->GetLast();j++){
1129 Double_t yval = yax->GetBinCenter(j);
1130 for(Int_t k=zax->GetFirst();k<=zax->GetLast();k++){
1131 Double_t zval = zax->GetBinCenter(k);
1132 Int_t bc =(Int_t)histo->GetBinContent(i,j,k);
1133 if(bc==0) continue;
1134 if(bc>nOfTracklets){
1135 nOfTracklets=bc;
1136 peak[2] = zval;
1137 peak[1] = yval;
1138 peak[0] = xval;
1139 peakbin[2] = k;
1140 peakbin[1] = j;
1141 peakbin[0] = i;
1142 peak2bin[2] = -1;
1143 peak2bin[1] = -1;
1144 peak2bin[0] = -1;
1145 bc2=-1;
1146 nOfTimes = 1;
1147 }else if(bc==nOfTracklets){
1148 if(TMath::Abs(i-peakbin[0])<=1 && TMath::Abs(j-peakbin[1])<=1 && TMath::Abs(k-peakbin[2])<=1){
1149 peak2bin[2] = k;
1150 peak2bin[1] = j;
1151 peak2bin[0] = i;
1152 bc2=bc;
1153 nOfTimes = 1;
1154 }else{
1155 nOfTimes++;
1156 }
1157 }
1158 }
1159 }
1160 }
1161 if(peak2bin[0]>=-1 && bc2!=-1){ // two contiguous peak-cells with same contents
1162 peak[0]=0.5*(xax->GetBinCenter(peakbin[0])+xax->GetBinCenter(peak2bin[0]));
1163 peak[1]=0.5*(yax->GetBinCenter(peakbin[1])+yax->GetBinCenter(peak2bin[1]));
1164 peak[2]=0.5*(zax->GetBinCenter(peakbin[2])+zax->GetBinCenter(peak2bin[2]));
1165 nOfTracklets+=bc2;
1166 nOfTimes=1;
1167 }
1168}
1169//________________________________________________________
1170void AliITSVertexer3D::MarkUsedClusters(){
1171 // Mark clusters of tracklets used in vertex claulation
1172 for(Int_t i=0; i<fLines.GetEntriesFast();i++){
1173 AliStrLine *lin = (AliStrLine*)fLines.At(i);
1174 Int_t idClu1=lin->GetIdPoint(0);
1175 Int_t idClu2=lin->GetIdPoint(1);
1176 fUsedCluster.SetBitNumber(idClu1);
1177 fUsedCluster.SetBitNumber(idClu2);
1178 }
1179}
1180//________________________________________________________
1181Int_t AliITSVertexer3D::RemoveTracklets(){
1182 // Remove trackelts close to first found vertex
1183 Double_t vert[3]={fVert3D.GetXv(),fVert3D.GetYv(),fVert3D.GetZv()};
1184 Int_t nRemoved=0;
1185 for(Int_t i=0; i<fLines.GetEntriesFast();i++){
1186 AliStrLine *lin = (AliStrLine*)fLines.At(i);
1187 if(lin->GetDistFromPoint(vert)<fDCAforPileup){
1188 Int_t idClu1=lin->GetIdPoint(0);
1189 Int_t idClu2=lin->GetIdPoint(1);
1190 fUsedCluster.SetBitNumber(idClu1);
1191 fUsedCluster.SetBitNumber(idClu2);
1192 fLines.RemoveAt(i);
1193 ++nRemoved;
1194 }
1195 }
1196 fLines.Compress();
1197 return nRemoved;
1198}
1199//________________________________________________________
1200void AliITSVertexer3D::FindOther3DVertices(TTree *itsClusterTree){
1201 // pileup identification based on 3D vertexing with not used clusters
1202
1203 fVertArray = new AliESDVertex[kMaxPileupVertices+1];
1204 fVertArray[0]=(*fCurrentVertex);
1205 Int_t nFoundVert=1;
1206 for(Int_t iPilV=1; iPilV<=kMaxPileupVertices; iPilV++){
1207 MarkUsedClusters();
1208 fLines.Clear("C");
1209 Int_t nolines = FindTracklets(itsClusterTree,2);
1210 if(nolines>=2){
1211 Int_t nr=RemoveTracklets();
1212 nolines-=nr;
1213 if(nolines>=2){
1214 Int_t rc=Prepare3DVertex(2);
1215 if(rc==0){
1216 fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
1217 if(fVert3D.GetNContributors()>=fMinTrackletsForPilup){
1218 fIsPileup=kTRUE;
1219 fVertArray[nFoundVert]=fVert3D;
1220 nFoundVert++;
1221 if(nFoundVert==2){
1222 fZpuv=fVert3D.GetZv();
1223 fNTrpuv=fVert3D.GetNContributors();
1224 }
1225 }
1226 }
1227 }
1228 }
1229 }
1230 fNoVertices=nFoundVert;
1231}
1232//______________________________________________________________________
1233void AliITSVertexer3D::PileupFromZ(){
1234 // Calls the pileup algorithm of ALiITSVertexerZ
1235 Int_t binmin, binmax;
1236 Int_t nPeaks=AliITSVertexerZ::GetPeakRegion(fZHisto,binmin,binmax);
1237 if(nPeaks==2)AliWarning("2 peaks found");
1238 Int_t firstPeakCont=0;
1239 Double_t firstPeakPos=0.;
1240 for(Int_t i=binmin-1;i<=binmax+1;i++){
1241 firstPeakCont+=static_cast<Int_t>(fZHisto->GetBinContent(i));
1242 firstPeakPos+=fZHisto->GetBinContent(i)*fZHisto->GetBinCenter(i);
1243 }
1244 if(firstPeakCont>0){
1245 firstPeakPos/=firstPeakCont;
1246 Int_t ncontr2=0;
1247 if(firstPeakCont>fMinTrackletsForPilup){
1248 Float_t secPeakPos;
1249 ncontr2=AliITSVertexerZ::FindSecondPeak(fZHisto,binmin,binmax,secPeakPos);
1250 if(ncontr2>=fMinTrackletsForPilup){
1251 fIsPileup=kTRUE;
1252 fNoVertices=2;
1253 AliESDVertex secondVert(secPeakPos,0.1,ncontr2);
1254 fVertArray = new AliESDVertex[2];
1255 fVertArray[0]=(*fCurrentVertex);
1256 fVertArray[1]=secondVert;
1257 fZpuv=secPeakPos;
1258 fNTrpuv=ncontr2;
1259 }
1260 }
1261 }
1262}
1263
1264//________________________________________________________
1265Double_t AliITSVertexer3D::GetFraction(Int_t itr) const {
1266 // this method is used to fill a 3D histogram representing
1267 // the trajectories of the candidate tracklets
1268 // The computed fraction is used as a weight at filling time
1269 AliStrLine *str = (AliStrLine*)fLines.At(itr);
1270 Double_t spigolo=10.;
1271 Double_t cd[3];
1272 str->GetCd(cd);
1273 Double_t par=0.;
1274 Double_t maxl=TMath::Sqrt(3.)*spigolo;
1275 // intersection with a plane normal to the X axis
1276 if(TMath::AreEqualAbs(cd[0],0.,1.e-9)){
1277 par=1000000.;
1278 }
1279 else {
1280 par=spigolo/cd[0];
1281 }
1282 Double_t zc=cd[2]*par;
1283 Double_t yc=cd[1]*par;
1284 if((-spigolo<=yc && yc<=spigolo) && (-spigolo<=zc && zc<=spigolo))return TMath::Abs(par/maxl);
1285 // intersection with a plane normal to the Y axis
1286 if(TMath::AreEqualAbs(cd[1],0.,1.e-9)){
1287 par=1000000.;
1288 }
1289 else {
1290 par=spigolo/cd[1];
1291 }
1292 zc=cd[2]*par;
1293 Double_t xc=cd[0]*par;
1294 if((-spigolo<=xc && xc<=spigolo) && (-spigolo<=zc && zc<=spigolo))return TMath::Abs(par/maxl);
1295 // intersection with a plane normal to the Z axis
1296 if(TMath::AreEqualAbs(cd[2],0.,1.e-9)){
1297 par=1000000.;
1298 }
1299 else {
1300 par=spigolo/cd[2];
1301 }
1302 yc=cd[1]*par;
1303 xc=cd[0]*par;
1304 if((-spigolo<=xc && xc<=spigolo) && (-spigolo<=yc && yc<=spigolo))return TMath::Abs(par/maxl);
1305 // control should never reach the following lines
1306 AliError(Form("anomalous tracklet direction for tracklet %d in fLines\n",itr));
1307 str->PrintStatus();
1308 return 0.;
1309}
1310
1311//________________________________________________________
1312void AliITSVertexer3D::PrintStatus() const {
1313 // Print current status
1314 printf("========= First step selections =====================\n");
1315 printf("Cut on diamond (Z) %f\n",fZCutDiamond);
1316 printf("Loose cut on Delta Phi %f\n",fCoarseDiffPhiCut);
1317 printf("Loose cut on tracklet DCA to Z axis %f\n",fCoarseMaxRCut);
1318 printf("Cut on DCA - tracklet to tracklet and to vertex %f\n",fDCAcut);
1319 printf("========= Second step selections ====================\n");
1320 printf("Cut on tracklet-to-first-vertex Z distance %f\n",fMaxZCut);
1321 printf("Max Phi difference: %f\n",fDiffPhiMax);
1322 printf("Cut on tracklet DCA to beam axis %f\n",fMaxRCut);
1323 printf("Cut on tracklet DCA to beam axis (algo2) %f\n",fMaxRCut2);
1324 printf("========= Pileup selections =========================\n");
1325 printf("Pileup algo: %d\n",fPileupAlgo);
1326 printf("Min DCA to 1st vertex for pileup (algo 0 and 1): %f\n",fDCAforPileup);
1327 printf("Cut on distance between pair-vertices (algo 2): %f\n",fCutOnPairs);
1328 printf("Maximum number of clusters on L1 or L2 for downscale: %d\n",fMaxNumOfClForDownScale);
1329 printf("Maximum number of clusters on L1 or L2 for histo rebin: %d\n",fMaxNumOfClForRebin);
1330 printf("=======================================================\n");
1331}
1332