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