ITS/AliITSVertexer3D.cxx

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