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