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