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