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