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 "AliESDVertex.h"
  17 #include "AliLog.h"
  18 #include "AliStrLine.h"
  19 #include "AliTracker.h"
  20 #include "AliITSDetTypeRec.h"
  21 #include "AliITSRecPoint.h"
  22 #include "AliITSgeomTGeo.h"
  23 #include "AliVertexerTracks.h"
  24 #include "AliITSVertexer3D.h"
  25 #include "AliITSVertexerZ.h"
  26 #include "AliITSSortTrkl.h"
  27 /////////////////////////////////////////////////////////////////
  28 // this class implements a method to determine
  29 // the 3 coordinates of the primary vertex
  30 // for p-p collisions
  31 // It can be used successfully with Pb-Pb collisions
  32 ////////////////////////////////////////////////////////////////
  33
  34 ClassImp(AliITSVertexer3D)
  35
  36 /* $Id$ */
  37
  38 //______________________________________________________________________
  39 AliITSVertexer3D::AliITSVertexer3D():
  40   AliITSVertexer(),
  41   fLines("AliStrLine",1000),
  42   fVert3D(),
  43   fCoarseDiffPhiCut(0.),
  44   fFineDiffPhiCut(0.),
  45   fCutOnPairs(0.),
  46   fCoarseMaxRCut(0.),
  47   fMaxRCut(0.),
  48   fZCutDiamond(0.),
  49   fMaxZCut(0.),
  50   fDCAcut(0.),
  51   fDiffPhiMax(0.),
  52   fMeanPSelTrk(0.),
  53   fMeanPtSelTrk(0.),
  54   fUsedCluster(kMaxCluPerMod*kNSPDMod),
  55   fZHisto(0),
  56   fDCAforPileup(0.),
  57   fDiffPhiforPileup(0.),
  58   fBinSizeR(0.),
  59   fBinSizeZ(0.),
  60   fPileupAlgo(0)
  61 {
  62   // Default constructor
  63   SetCoarseDiffPhiCut();
  64   SetFineDiffPhiCut();
  65   SetCutOnPairs();
  66   SetCoarseMaxRCut();
  67   SetMaxRCut();
  68   SetZCutDiamond();
  69   SetMaxZCut();
  70   SetDCACut();
  71   SetDiffPhiMax();
  72   SetMeanPSelTracks();
  73   SetMeanPtSelTracks();
  74   SetMinDCAforPileup();
  75   SetDeltaPhiforPileup();
  76   SetPileupAlgo();
  77   SetBinSizeR();
  78   SetBinSizeZ();
  79   Double_t binsize=0.02; // default 200 micron
  80   Int_t nbins=static_cast<Int_t>(1+2*fZCutDiamond/binsize);
  81   fZHisto=new TH1F("hz","",nbins,-fZCutDiamond,-fZCutDiamond+binsize*nbins);
  82 }
  83
  84 //______________________________________________________________________
  85 AliITSVertexer3D::~AliITSVertexer3D() {
  86   // Destructor
  87   fLines.Clear("C");
  88   if(fZHisto) delete fZHisto;
  89 }
  90
  91 //______________________________________________________________________
  92 void AliITSVertexer3D::ResetVert3D(){
  93   //
  94   ResetVertex();
  95   fVert3D.SetXv(0.);
  96   fVert3D.SetYv(0.);
  97   fVert3D.SetZv(0.);
  98   fVert3D.SetDispersion(0.);
  99   fVert3D.SetNContributors(0);
 100   fUsedCluster.ResetAllBits(0);
 101 }
 102 //______________________________________________________________________
 103 AliESDVertex* AliITSVertexer3D::FindVertexForCurrentEvent(TTree *itsClusterTree){
 104   // Defines the AliESDVertex for the current event
 105   ResetVert3D();
 106   AliDebug(1,"FindVertexForCurrentEvent - 3D - PROCESSING NEXT EVENT");
 107   fLines.Clear("C");
 108   fCurrentVertex = NULL;
 109
 110   Int_t nolines = FindTracklets(itsClusterTree,0);
 111   if(nolines>=2){
 112     Int_t rc=Prepare3DVertex(0);
 113     if(fPileupAlgo == 2 && rc == 0) FindVertex3DIterative();
 114     else if(fPileupAlgo<2 && rc == 0) FindVertex3D(itsClusterTree);
 115   }
 116
 117   if(!fCurrentVertex){
 118     AliITSVertexerZ vertz(GetNominalPos()[0],GetNominalPos()[1]);
 119     vertz.SetDetTypeRec(GetDetTypeRec());
 120     AliDebug(1,"Call Vertexer Z\n");
 121     vertz.SetLowLimit(-fZCutDiamond);
 122     vertz.SetHighLimit(fZCutDiamond);
 123     AliESDVertex* vtxz = vertz.FindVertexForCurrentEvent(itsClusterTree);
 124     if(vtxz){
 125       Double_t position[3]={GetNominalPos()[0],GetNominalPos()[1],vtxz->GetZv()};
 126       Double_t covmatrix[6];
 127       vtxz->GetCovMatrix(covmatrix);
 128       Double_t chi2=99999.;
 129       Int_t    nContr=vtxz->GetNContributors();
 130       fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr);
 131       fCurrentVertex->SetTitle("vertexer: Z");
 132       fCurrentVertex->SetName("SPDVertexZ");
 133       delete vtxz;
 134     }
 135
 136   }
 137   FindMultiplicity(itsClusterTree);
 138   return fCurrentVertex;
 139 }
 140
 141 //______________________________________________________________________
 142 void AliITSVertexer3D::FindVertex3D(TTree *itsClusterTree){
 143   // 3D algorithm
 144   /*  uncomment to debug
 145       printf("Vertex found in first iteration:\n");
 146       fVert3D.Print();
 147       printf("Start second iteration\n");
 148       end of debug lines  */
 149   if(fVert3D.GetNContributors()>0){
 150     fLines.Clear("C");
 151     Int_t nolines = FindTracklets(itsClusterTree,1);
 152     if(nolines>=2){
 153       Int_t rc=Prepare3DVertex(1);
 154       if(rc!=0) fVert3D.SetNContributors(0); // exclude this vertex
 155     }
 156   }
 157   /*  uncomment to debug
 158       printf("Vertex found in second iteration:\n");
 159       fVert3D.Print();
 160       end of debug lines  */
 161
 162   Double_t vRadius=TMath::Sqrt(fVert3D.GetXv()*fVert3D.GetXv()+fVert3D.GetYv()*fVert3D.GetYv());
 163   if(vRadius<GetPipeRadius() && fVert3D.GetNContributors()>0){
 164     Double_t position[3]={fVert3D.GetXv(),fVert3D.GetYv(),fVert3D.GetZv()};
 165     Double_t covmatrix[6];
 166     fVert3D.GetCovMatrix(covmatrix);
 167     Double_t chi2=99999.;
 168     Int_t    nContr=fVert3D.GetNContributors();
 169     fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr);
 170     fCurrentVertex->SetTitle("vertexer: 3D");
 171     fCurrentVertex->SetName("SPDVertex3D");
 172     fCurrentVertex->SetDispersion(fVert3D.GetDispersion());
 173     fNoVertices=1;
 174
 175     switch(fPileupAlgo){
 176     case 0: PileupFromZ(); break;
 177     case 1: FindOther3DVertices(itsClusterTree); break;
 178     default: AliError("Wrong pileup algorithm"); break;
 179     }
 180     if(fNoVertices==1){
 181       fVertArray = new AliESDVertex[1];
 182       fVertArray[0]=(*fCurrentVertex);
 183     }
 184   }
 185 }
 186
 187 //______________________________________________________________________
 188 void AliITSVertexer3D::FindVertex3DIterative(){
 189   // Defines the AliESDVertex for the current event
 190   Int_t numsor=fLines.GetEntriesFast()*(fLines.GetEntriesFast()-1)/2;
 191   //cout<<"AliITSVertexer3D::FindVertexForCurentEvent: Number of tracklets selected for vertexing "<<fLines.GetEntriesFast()<<"; Number of pairs: "<<numsor<<endl;
 192   AliITSSortTrkl srt(fLines,numsor,fCutOnPairs,fCoarseMaxRCut);
 193   srt.FindClusters();
 194   AliInfo(Form("Number of vertices: %d",srt.GetNumberOfClusters()));
 195
 196   fNoVertices = srt.GetNumberOfClusters();
 197   //printf("fNoVertices = %d \n",fNoVertices);
 198   if(fNoVertices>0){
 199     fVertArray = new AliESDVertex[fNoVertices];
 200     for(Int_t kk=0; kk<srt.GetNumberOfClusters(); kk++){
 201       Int_t size = 0;
 202       Int_t *labels = srt.GetTrackletsLab(kk,size);
 203       /*
 204         Int_t *pairs = srt.GetClusters(kk);
 205         Int_t nopai = srt.GetSizeOfCluster(kk);
 206         cout<<"***** Vertex number "<<kk<<".  Pairs: \n";
 207         for(Int_t jj=0;jj<nopai;jj++){
 208         cout<<pairs[jj]<<" - ";
 209         if(jj>0 & jj%8==0)cout<<endl;
 210         }
 211         cout<<endl;
 212         cout<<"***** Vertex number "<<kk<<".  Labels: \n";
 213       */
 214       AliStrLine **tclo = new AliStrLine* [size];
 215       for(Int_t jj=0;jj<size;jj++){
 216         //          cout<<labels[jj]<<" - ";
 217         //          if(jj>0 & jj%8==0)cout<<endl;
 218         tclo[jj] = dynamic_cast<AliStrLine*>(fLines[labels[jj]]);
 219       }
 220       //          cout<<endl;
 221       delete []labels;
 222       fVertArray[kk]=AliVertexerTracks::TrackletVertexFinder(tclo,size);
 223       delete [] tclo;
 224       //          fVertArray[kk].PrintStatus();
 225       if(kk == 1){
 226         // at least one second vertex is present
 227         fIsPileup = kTRUE;
 228         fNTrpuv = fVertArray[kk].GetNContributors();
 229         fZpuv = fVertArray[kk].GetZv();
 230       }
 231     }
 232     Double_t vRadius=TMath::Sqrt(fVertArray[0].GetXv()*fVertArray[0].GetXv()+fVertArray[0].GetYv()*fVertArray[0].GetYv());
 233     if(vRadius<GetPipeRadius() && fVertArray[0].GetNContributors()>0){
 234       Double_t position[3]={fVertArray[0].GetXv(),fVertArray[0].GetYv(),fVertArray[0].GetZv()};
 235       Double_t covmatrix[6];
 236       fVertArray[0].GetCovMatrix(covmatrix);
 237       Double_t chi2=99999.;
 238       Int_t    nContr=fVertArray[0].GetNContributors();
 239       fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr);
 240       fCurrentVertex->SetTitle("vertexer: 3D");
 241       fCurrentVertex->SetName("SPDVertex3D");
 242       fCurrentVertex->SetDispersion(fVertArray[0].GetDispersion());
 243     }
 244   }
 245
 246 }
 247
 248 //______________________________________________________________________
 249 Bool_t AliITSVertexer3D::DistBetweenVertices(AliESDVertex &a, AliESDVertex &b, Double_t test, Double_t &dist){
 250   // method to compare the distance between vertices a and b with "test"
 251   //it returns kTRUE is the distance is less or equal to test
 252   dist = (a.GetX()-b.GetX()) * (a.GetX()-b.GetX());
 253   dist +=  (a.GetY()-b.GetY()) * (a.GetY()-b.GetY());
 254   dist +=  (a.GetZ()-b.GetZ()) * (a.GetZ()-b.GetZ());
 255   dist = TMath::Sqrt(dist);
 256   if(dist <= test)return kTRUE;
 257   return kFALSE;
 258 }
 259
 260
 261 //______________________________________________________________________
 262 Int_t AliITSVertexer3D::FindTracklets(TTree *itsClusterTree, Int_t optCuts){
 263   // All the possible combinations between recpoints on layer 1and 2 are
 264   // considered. Straight lines (=tracklets)are formed.
 265   // The tracklets are processed in Prepare3DVertex
 266
 267   if(!GetDetTypeRec())AliFatal("DetTypeRec pointer has not been set");
 268
 269   TTree *tR = itsClusterTree;
 270   fDetTypeRec->ResetRecPoints();
 271   fDetTypeRec->SetTreeAddressR(tR);
 272   TClonesArray *itsRec  = 0;
 273   if(optCuts==0) fZHisto->Reset();
 274  // gc1 are local and global coordinates for layer 1
 275   Float_t gc1f[3]={0.,0.,0.};
 276   Double_t gc1[3]={0.,0.,0.};
 277   // gc2 are local and global coordinates for layer 2
 278   Float_t gc2f[3]={0.,0.,0.};
 279   Double_t gc2[3]={0.,0.,0.};
 280
 281   itsRec = fDetTypeRec->RecPoints();
 282   TBranch *branch = NULL;
 283   branch = tR->GetBranch("ITSRecPoints");
 284   if(!branch){
 285     AliWarning("Null pointer for RecPoints branch, vertex not calculated");
 286     return -1;
 287   }
 288
 289   // Set values for cuts
 290   Double_t xbeam=GetNominalPos()[0];
 291   Double_t ybeam=GetNominalPos()[1];
 292   Double_t zvert=0.;
 293   Double_t deltaPhi=fCoarseDiffPhiCut;
 294   Double_t deltaR=fCoarseMaxRCut;
 295   Double_t dZmax=fZCutDiamond;
 296   if(fPileupAlgo == 2){
 297     deltaPhi=fFineDiffPhiCut;
 298     deltaR=fMaxRCut;
 299     if(optCuts != 0)AliWarning(Form("fPileupAlgo=2 AND optCuts=%d has been selected. It should be 0",optCuts));
 300   } else if(optCuts==1){
 301     xbeam=fVert3D.GetXv();
 302     ybeam=fVert3D.GetYv();
 303     zvert=fVert3D.GetZv();
 304     deltaPhi = fDiffPhiMax;
 305     deltaR=fMaxRCut;
 306     dZmax=fMaxZCut;
 307   } else if(optCuts==2){
 308     xbeam=fVert3D.GetXv();
 309     ybeam=fVert3D.GetYv();
 310     deltaPhi = fDiffPhiforPileup;
 311     deltaR=fMaxRCut;
 312   }
 313
 314   Int_t nrpL1 = 0;    // number of rec points on layer 1
 315   Int_t nrpL2 = 0;    // number of rec points on layer 2
 316
 317   // By default irstL1=0 and lastL1=79
 318   Int_t firstL1 = AliITSgeomTGeo::GetModuleIndex(1,1,1);
 319   Int_t lastL1 = AliITSgeomTGeo::GetModuleIndex(2,1,1)-1;
 320   for(Int_t module= firstL1; module<=lastL1;module++){  // count number of recopints on layer 1
 321     branch->GetEvent(module);
 322     nrpL1+= itsRec->GetEntries();
 323     fDetTypeRec->ResetRecPoints();
 324   }
 325   //By default firstL2=80 and lastL2=239
 326   Int_t firstL2 = AliITSgeomTGeo::GetModuleIndex(2,1,1);
 327   Int_t lastL2 = AliITSgeomTGeo::GetModuleIndex(3,1,1)-1;
 328   for(Int_t module= firstL2; module<=lastL2;module++){  // count number of recopints on layer 2
 329     branch->GetEvent(module);
 330     nrpL2+= itsRec->GetEntries();
 331     fDetTypeRec->ResetRecPoints();
 332   }
 333   if(nrpL1 == 0 || nrpL2 == 0){
 334     AliDebug(1,Form("No RecPoints in at least one SPD layer (%d %d)",nrpL1,nrpL2));
 335     return -1;
 336   }
 337   AliDebug(1,Form("RecPoints on Layer 1,2 = %d, %d\n",nrpL1,nrpL2));
 338
 339   Double_t a[3]={xbeam,ybeam,0.};
 340   Double_t b[3]={xbeam,ybeam,10.};
 341   AliStrLine zeta(a,b,kTRUE);
 342   static Double_t bField=TMath::Abs(AliTracker::GetBz()/10.); //T
 343   SetMeanPPtSelTracks(bField);
 344
 345   Int_t nolines = 0;
 346   // Loop on modules of layer 1
 347   for(Int_t modul1= firstL1; modul1<=lastL1;modul1++){   // Loop on modules of layer 1
 348     if(!fUseModule[modul1]) continue;
 349     UShort_t ladder=int(modul1/4)+1; // ladders are numbered starting from 1
 350     branch->GetEvent(modul1);
 351     Int_t nrecp1 = itsRec->GetEntries();
 352     static TClonesArray prpl1("AliITSRecPoint",nrecp1);
 353     prpl1.SetOwner();
 354     for(Int_t j=0;j<nrecp1;j++){
 355       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
 356       new(prpl1[j])AliITSRecPoint(*recp);
 357     }
 358     fDetTypeRec->ResetRecPoints();
 359     for(Int_t j=0;j<nrecp1;j++){
 360       if(j>kMaxCluPerMod) continue;
 361       UShort_t idClu1=modul1*kMaxCluPerMod+j;
 362       if(fUsedCluster.TestBitNumber(idClu1)) continue;
 363       AliITSRecPoint *recp1 = (AliITSRecPoint*)prpl1.At(j);
 364       recp1->GetGlobalXYZ(gc1f);
 365       for(Int_t ico=0;ico<3;ico++)gc1[ico]=gc1f[ico];
 366
 367       Double_t phi1 = TMath::ATan2(gc1[1]-ybeam,gc1[0]-xbeam);
 368       if(phi1<0)phi1=2*TMath::Pi()+phi1;
 369       for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){
 370         for(Int_t k=0;k<4;k++){
 371           Int_t ladmod=fLadders[ladder-1]+ladl2;
 372           if(ladmod>AliITSgeomTGeo::GetNLadders(2)) ladmod=ladmod-AliITSgeomTGeo::GetNLadders(2);
 373           Int_t modul2=AliITSgeomTGeo::GetModuleIndex(2,ladmod,k+1);
 374           if(!fUseModule[modul2]) continue;
 375           branch->GetEvent(modul2);
 376           Int_t nrecp2 = itsRec->GetEntries();
 377           for(Int_t j2=0;j2<nrecp2;j2++){
 378             if(j2>kMaxCluPerMod) continue;
 379             UShort_t idClu2=modul2*kMaxCluPerMod+j2;
 380             if(fUsedCluster.TestBitNumber(idClu2)) continue;
 381             AliITSRecPoint *recp2 = (AliITSRecPoint*)itsRec->At(j2);
 382             recp2->GetGlobalXYZ(gc2f);
 383             for(Int_t ico=0;ico<3;ico++)gc2[ico]=gc2f[ico];
 384             Double_t phi2 = TMath::ATan2(gc2[1]-ybeam,gc2[0]-xbeam);
 385             if(phi2<0)phi2=2*TMath::Pi()+phi2;
 386             Double_t diff = TMath::Abs(phi2-phi1);
 387             if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
 388             if(optCuts==0 && diff<fDiffPhiMax){
 389               Double_t r1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]);
 390               Double_t zc1=gc1[2];
 391               Double_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
 392               Double_t zc2=gc2[2];
 393               Double_t zr0=(r2*zc1-r1*zc2)/(r2-r1); //Z @ null radius
 394               fZHisto->Fill(zr0);
 395             }
 396             if(diff>deltaPhi)continue;
 397             AliStrLine line(gc1,gc2,kTRUE);
 398             Double_t cp[3];
 399             Int_t retcode = line.Cross(&zeta,cp);
 400             if(retcode<0)continue;
 401             Double_t dca = line.GetDCA(&zeta);
 402             if(dca<0.) continue;
 403             if(dca>deltaR)continue;
 404             Double_t deltaZ=cp[2]-zvert;
 405             if(TMath::Abs(deltaZ)>dZmax)continue;
 406
 407             if(nolines == 0){
 408               if(fLines.GetEntriesFast()>0)fLines.Clear("C");
 409             }
 410             Float_t cov[6];
 411             recp2->GetGlobalCov(cov);
 412
 413
 414             Double_t rad1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]);
 415             Double_t rad2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
 416             Double_t factor=(rad1+rad2)/(rad2-rad1); //factor to account for error on tracklet direction
 417
 418             Double_t curvErr=0;
 419             if(bField>0.00001){
 420               Double_t curvRadius=fMeanPtSelTrk/(0.3*bField)*100; //cm
 421               Double_t dRad=TMath::Sqrt((gc1[0]-gc2[0])*(gc1[0]-gc2[0])+(gc1[1]-gc2[1])*(gc1[1]-gc2[1]));
 422               Double_t aux=dRad/2.+rad1;
 423               curvErr=TMath::Sqrt(curvRadius*curvRadius-dRad*dRad/4.)-TMath::Sqrt(curvRadius*curvRadius-aux*aux); //cm
 424             }
 425             Double_t sigmasq[3];
 426             sigmasq[0]=(cov[0]+curvErr*curvErr/2.)*factor*factor;
 427             sigmasq[1]=(cov[3]+curvErr*curvErr/2.)*factor*factor;
 428             sigmasq[2]=cov[5]*factor*factor;
 429
 430             // Multiple scattering
 431             Double_t pOverMass=fMeanPSelTrk/0.140;
 432             Double_t beta2=pOverMass*pOverMass/(1+pOverMass*pOverMass);
 433             Double_t p2=fMeanPSelTrk*fMeanPSelTrk;
 434             Double_t rBP=GetPipeRadius();
 435             Double_t dBP=0.08/35.3; // 800 um of Be
 436             Double_t dL1=0.01; //approx. 1% of radiation length
 437             Double_t theta2BP=14.1*14.1/(beta2*p2*1e6)*dBP;
 438             Double_t theta2L1=14.1*14.1/(beta2*p2*1e6)*dL1;
 439             Double_t rtantheta1=(rad2-rad1)*TMath::Tan(TMath::Sqrt(theta2L1));
 440             Double_t rtanthetaBP=(rad1-rBP)*TMath::Tan(TMath::Sqrt(theta2BP));
 441             for(Int_t ico=0; ico<3;ico++){
 442               sigmasq[ico]+=rtantheta1*rtantheta1*factor*factor/3.;
 443               sigmasq[ico]+=rtanthetaBP*rtanthetaBP*factor*factor/3.;
 444             }
 445             Double_t wmat[9]={1.,0.,0.,0.,1.,0.,0.,0.,1.};
 446             if(sigmasq[0]!=0.) wmat[0]=1./sigmasq[0];
 447             if(sigmasq[1]!=0.) wmat[4]=1./sigmasq[1];
 448             if(sigmasq[2]!=0.) wmat[8]=1./sigmasq[2];
 449             new(fLines[nolines++])AliStrLine(gc1,sigmasq,wmat,gc2,kTRUE,idClu1,idClu2);
 450
 451           }
 452           fDetTypeRec->ResetRecPoints();
 453         }
 454       }
 455     }
 456     prpl1.Clear();
 457   }
 458   if(nolines == 0)return -2;
 459   return nolines;
 460 }
 461
 462 //______________________________________________________________________
 463 Int_t  AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){
 464   // Finds the 3D vertex information using tracklets
 465   Int_t retcode = -1;
 466
 467   Double_t xbeam=GetNominalPos()[0];
 468   Double_t ybeam=GetNominalPos()[1];
 469   Double_t zvert=0.;
 470   Double_t deltaR=fCoarseMaxRCut;
 471   if(fPileupAlgo == 2) {
 472     deltaR=fMaxRCut;
 473     if(optCuts!=0)AliWarning(Form("fPileupAlgo=2 AND optCuts=%d. It should be 0",optCuts));
 474   }
 475   Double_t dZmax=fZCutDiamond;
 476   if(optCuts==1){
 477     xbeam=fVert3D.GetXv();
 478     ybeam=fVert3D.GetYv();
 479     zvert=fVert3D.GetZv();
 480     deltaR=fMaxRCut;
 481     dZmax=fMaxZCut;
 482   }else if(optCuts==2){
 483     xbeam=fVert3D.GetXv();
 484     ybeam=fVert3D.GetYv();
 485     deltaR=fMaxRCut;
 486   }
 487
 488   Double_t rl=-fCoarseMaxRCut;
 489   Double_t rh=fCoarseMaxRCut;
 490   Double_t zl=-fZCutDiamond;
 491   Double_t zh=fZCutDiamond;
 492   Int_t nbr=(Int_t)((rh-rl)/fBinSizeR+0.0001);
 493   Int_t nbz=(Int_t)((zh-zl)/fBinSizeZ+0.0001);
 494   Int_t nbrcs=(Int_t)((rh-rl)/(fBinSizeR*2.)+0.0001);
 495   Int_t nbzcs=(Int_t)((zh-zl)/(fBinSizeZ*2.)+0.0001);
 496
 497   TH3F *h3d = NULL;
 498   TH3F *h3dcs = NULL;
 499   if(fPileupAlgo !=2){
 500     h3d = new TH3F("h3d","xyz distribution",nbr,rl,rh,nbr,rl,rh,nbz,zl,zh);
 501     h3dcs = new TH3F("h3dcs","xyz distribution",nbrcs,rl,rh,nbrcs,rl,rh,nbzcs,zl,zh);
 502   }
 503   // cleanup of the TCLonesArray of tracklets (i.e. fakes are removed)
 504   Int_t *validate = new Int_t [fLines.GetEntriesFast()];
 505   for(Int_t i=0; i<fLines.GetEntriesFast();i++)validate[i]=0;
 506   for(Int_t i=0; i<fLines.GetEntriesFast()-1;i++){
 507     AliStrLine *l1 = (AliStrLine*)fLines.At(i);
 508     for(Int_t j=i+1;j<fLines.GetEntriesFast();j++){
 509       AliStrLine *l2 = (AliStrLine*)fLines.At(j);
 510       Double_t dca=l1->GetDCA(l2);
 511       if(dca > fDCAcut || dca<0.00001) continue;
 512       Double_t point[3];
 513       Int_t retc = l1->Cross(l2,point);
 514       if(retc<0)continue;
 515       Double_t deltaZ=point[2]-zvert;
 516      if(TMath::Abs(deltaZ)>dZmax)continue;
 517       Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
 518       if(rad>fCoarseMaxRCut)continue;
 519       Double_t deltaX=point[0]-xbeam;
 520       Double_t deltaY=point[1]-ybeam;
 521       Double_t raddist=TMath::Sqrt(deltaX*deltaX+deltaY*deltaY);
 522       if(raddist>deltaR)continue;
 523       validate[i]=1;
 524       validate[j]=1;
 525       if(fPileupAlgo != 2){
 526         h3d->Fill(point[0],point[1],point[2]);
 527         h3dcs->Fill(point[0],point[1],point[2]);
 528       }
 529     }
 530   }
 531
 532
 533
 534   Int_t numbtracklets=0;
 535   for(Int_t i=0; i<fLines.GetEntriesFast();i++)if(validate[i]>=1)numbtracklets++;
 536   if(numbtracklets<2){
 537     delete [] validate;
 538     if(fPileupAlgo != 2){
 539       delete h3d;
 540       delete h3dcs;
 541     }
 542     return retcode;
 543   }
 544
 545   for(Int_t i=0; i<fLines.GetEntriesFast();i++){
 546     if(validate[i]<1)fLines.RemoveAt(i);
 547   }
 548   fLines.Compress();
 549   AliDebug(1,Form("Number of tracklets (after compress)%d ",fLines.GetEntriesFast()));
 550   delete [] validate;
 551
 552   // Exit here if Pileup Algorithm 2 has been chosen
 553   if(fPileupAlgo == 2)return 0;
 554   //
 555
 556
 557   //        Find peaks in histos
 558
 559   Double_t peak[3]={0.,0.,0.};
 560   Int_t ntrkl,ntimes;
 561   FindPeaks(h3d,peak,ntrkl,ntimes);
 562   delete h3d;
 563   Double_t binsizer=(rh-rl)/nbr;
 564   Double_t binsizez=(zh-zl)/nbz;
 565   if(optCuts==0 && (ntrkl<=2 || ntimes>1)){
 566     ntrkl=0;
 567     ntimes=0;
 568     FindPeaks(h3dcs,peak,ntrkl,ntimes);
 569     binsizer=(rh-rl)/nbrcs;
 570     binsizez=(zh-zl)/nbzcs;
 571     if(ntrkl==1 || ntimes>1){delete h3dcs; return retcode;}
 572   }
 573   delete h3dcs;
 574
 575   //         Second selection loop
 576
 577   Double_t bs=(binsizer+binsizez)/2.;
 578   for(Int_t i=0; i<fLines.GetEntriesFast();i++){
 579     AliStrLine *l1 = (AliStrLine*)fLines.At(i);
 580     if(l1->GetDistFromPoint(peak)>2.5*bs)fLines.RemoveAt(i);
 581   }
 582   fLines.Compress();
 583   AliDebug(1,Form("Number of tracklets (after 2nd compression) %d",fLines.GetEntriesFast()));
 584
 585   if(fLines.GetEntriesFast()>1){
 586     retcode=0;
 587     //  find a first candidate for the primary vertex
 588     fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
 589     // make a further selection on tracklets based on this first candidate
 590     fVert3D.GetXYZ(peak);
 591     AliDebug(1,Form("FIRST V candidate: %f ; %f ; %f",peak[0],peak[1],peak[2]));
 592     for(Int_t i=0; i<fLines.GetEntriesFast();i++){
 593       AliStrLine *l1 = (AliStrLine*)fLines.At(i);
 594       if(l1->GetDistFromPoint(peak)> fDCAcut)fLines.RemoveAt(i);
 595     }
 596     fLines.Compress();
 597     AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines.GetEntriesFast()));
 598     if(fLines.GetEntriesFast()>1){// this new tracklet selection is used
 599       fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
 600     }
 601   }
 602   return retcode;
 603 }
 604
 605 //________________________________________________________
 606 void AliITSVertexer3D::SetMeanPPtSelTracks(Double_t fieldTesla){
 607   // Sets mean values of Pt based on the field
 608   // for P (used in multiple scattering) the most probable value is used
 609   if(TMath::Abs(fieldTesla-0.5)<0.01){
 610     SetMeanPSelTracks(0.375);
 611     SetMeanPtSelTracks(0.630);
 612   }else if(TMath::Abs(fieldTesla-0.4)<0.01){
 613     SetMeanPSelTracks(0.375);
 614     SetMeanPtSelTracks(0.580);
 615   }else if(TMath::Abs(fieldTesla-0.2)<0.01){
 616     SetMeanPSelTracks(0.375);
 617     SetMeanPtSelTracks(0.530);
 618   }else if(fieldTesla<0.00001){
 619     SetMeanPSelTracks(0.375);
 620     SetMeanPtSelTracks(0.230);
 621   }else{
 622     SetMeanPSelTracks();
 623     SetMeanPtSelTracks();
 624   }
 625 }
 626
 627 //________________________________________________________
 628 void AliITSVertexer3D::FindPeaks(TH3F* histo, Double_t *peak, Int_t &nOfTracklets, Int_t &nOfTimes){
 629   // Finds bin with max contents in 3D histo of tracket intersections
 630   TAxis *xax = histo->GetXaxis();
 631   TAxis *yax = histo->GetYaxis();
 632   TAxis *zax = histo->GetZaxis();
 633   peak[0]=0.;
 634   peak[1]=0.;
 635   peak[2]=0.;
 636   nOfTracklets = 0;
 637   nOfTimes=0;
 638   Int_t peakbin[3]={0,0,0};
 639   Int_t peak2bin[3]={-1,-1,-1};
 640   Int_t bc2=-1;
 641   for(Int_t i=xax->GetFirst();i<=xax->GetLast();i++){
 642     Double_t xval = xax->GetBinCenter(i);
 643     for(Int_t j=yax->GetFirst();j<=yax->GetLast();j++){
 644       Double_t yval = yax->GetBinCenter(j);
 645       for(Int_t k=zax->GetFirst();k<=zax->GetLast();k++){
 646         Double_t zval = zax->GetBinCenter(k);
 647         Int_t bc =(Int_t)histo->GetBinContent(i,j,k);
 648         if(bc==0) continue;
 649         if(bc>nOfTracklets){
 650           nOfTracklets=bc;
 651           peak[2] = zval;
 652           peak[1] = yval;
 653           peak[0] = xval;
 654           peakbin[2] = k;
 655           peakbin[1] = j;
 656           peakbin[0] = i;
 657           peak2bin[2] = -1;
 658           peak2bin[1] = -1;
 659           peak2bin[0] = -1;
 660           bc2=-1;
 661           nOfTimes = 1;
 662         }else if(bc==nOfTracklets){
 663           if(TMath::Abs(i-peakbin[0])<=1 && TMath::Abs(j-peakbin[1])<=1 && TMath::Abs(k-peakbin[2])<=1){
 664             peak2bin[2] = k;
 665             peak2bin[1] = j;
 666             peak2bin[0] = i;
 667             bc2=bc;
 668             nOfTimes = 1;
 669           }else{
 670             nOfTimes++;
 671           }
 672         }
 673       }
 674     }
 675   }
 676   if(peak2bin[0]>=-1 && bc2!=-1){ // two contiguous peak-cells with same contents
 677     peak[0]=0.5*(xax->GetBinCenter(peakbin[0])+xax->GetBinCenter(peak2bin[0]));
 678     peak[1]=0.5*(yax->GetBinCenter(peakbin[1])+yax->GetBinCenter(peak2bin[1]));
 679     peak[2]=0.5*(zax->GetBinCenter(peakbin[2])+zax->GetBinCenter(peak2bin[2]));
 680     nOfTracklets+=bc2;
 681     nOfTimes=1;
 682   }
 683 }
 684 //________________________________________________________
 685 void AliITSVertexer3D::MarkUsedClusters(){
 686   // Mark clusters of tracklets used in vertex claulation
 687   for(Int_t i=0; i<fLines.GetEntriesFast();i++){
 688     AliStrLine *lin = (AliStrLine*)fLines.At(i);
 689     Int_t idClu1=lin->GetIdPoint(0);
 690     Int_t idClu2=lin->GetIdPoint(1);
 691     fUsedCluster.SetBitNumber(idClu1);
 692     fUsedCluster.SetBitNumber(idClu2);
 693   }
 694 }
 695 //________________________________________________________
 696 Int_t AliITSVertexer3D::RemoveTracklets(){
 697   // Remove trackelts close to first found vertex
 698   Double_t vert[3]={fVert3D.GetXv(),fVert3D.GetYv(),fVert3D.GetZv()};
 699   Int_t nRemoved=0;
 700   for(Int_t i=0; i<fLines.GetEntriesFast();i++){
 701     AliStrLine *lin = (AliStrLine*)fLines.At(i);
 702     if(lin->GetDistFromPoint(vert)<fDCAforPileup){
 703       Int_t idClu1=lin->GetIdPoint(0);
 704       Int_t idClu2=lin->GetIdPoint(1);
 705       fUsedCluster.SetBitNumber(idClu1);
 706       fUsedCluster.SetBitNumber(idClu2);
 707       fLines.RemoveAt(i);
 708       ++nRemoved;
 709     }
 710   }
 711   fLines.Compress();
 712   return nRemoved;
 713 }
 714 //________________________________________________________
 715 void AliITSVertexer3D::FindOther3DVertices(TTree *itsClusterTree){
 716   // pileup identification based on 3D vertexing with not used clusters
 717   MarkUsedClusters();
 718   fLines.Clear("C");
 719   Int_t nolines = FindTracklets(itsClusterTree,2);
 720   if(nolines>=2){
 721     Int_t nr=RemoveTracklets();
 722     nolines-=nr;
 723     if(nolines>=2){
 724       Int_t rc=Prepare3DVertex(2);
 725       if(rc==0){
 726         fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
 727         if(fVert3D.GetNContributors()>=fMinTrackletsForPilup){
 728           fIsPileup=kTRUE;
 729           fNoVertices=2;
 730           fVertArray = new AliESDVertex[2];
 731           fVertArray[0]=(*fCurrentVertex);
 732           fVertArray[1]=fVert3D;
 733           fZpuv=fVert3D.GetZv();
 734           fNTrpuv=fVert3D.GetNContributors();
 735         }
 736       }
 737     }
 738   }
 739 }
 740 //______________________________________________________________________
 741 void AliITSVertexer3D::PileupFromZ(){
 742   // Calls the pileup algorithm of ALiITSVertexerZ
 743   Int_t binmin, binmax;
 744   Int_t nPeaks=AliITSVertexerZ::GetPeakRegion(fZHisto,binmin,binmax);
 745   if(nPeaks==2)AliWarning("2 peaks found");
 746   Int_t firstPeakCont=0;
 747   Double_t firstPeakPos=0.;
 748   for(Int_t i=binmin-1;i<=binmax+1;i++){
 749     firstPeakCont+=static_cast<Int_t>(fZHisto->GetBinContent(i));
 750     firstPeakPos+=fZHisto->GetBinContent(i)*fZHisto->GetBinCenter(i);
 751   }
 752   if(firstPeakCont>0){
 753     firstPeakPos/=firstPeakCont;
 754     Int_t ncontr2=0;
 755     if(firstPeakCont>fMinTrackletsForPilup){
 756       Float_t secPeakPos;
 757       ncontr2=AliITSVertexerZ::FindSecondPeak(fZHisto,binmin,binmax,secPeakPos);
 758       if(ncontr2>=fMinTrackletsForPilup){
 759         fIsPileup=kTRUE;
 760         fNoVertices=2;
 761         AliESDVertex secondVert(secPeakPos,0.1,ncontr2);
 762         fVertArray = new AliESDVertex[2];
 763         fVertArray[0]=(*fCurrentVertex);
 764         fVertArray[1]=secondVert;
 765         fZpuv=secPeakPos;
 766         fNTrpuv=ncontr2;
 767       }
 768     }
 769   }
 770 }
 771 //________________________________________________________
 772 void AliITSVertexer3D::PrintStatus() const {
 773   // Print current status
 774   printf("=======================================================\n");
 775   printf("Loose cut on Delta Phi %f\n",fCoarseDiffPhiCut);
 776   printf("Cut on tracklet DCA to Z axis %f\n",fCoarseMaxRCut);
 777   printf("Cut on tracklet DCA to beam axis %f\n",fMaxRCut);
 778   printf("Cut on diamond (Z) %f\n",fZCutDiamond);
 779   printf("Cut on DCA - tracklet to tracklet and to vertex %f\n",fDCAcut);
 780   printf("Max Phi difference: %f\n",fDiffPhiMax);
 781   printf("Pileup algo: %d\n",fPileupAlgo);
 782   printf("Min DCA to 1st vetrtex for pileup: %f\n",fDCAforPileup);
 783   printf("=======================================================\n");
 784 }