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