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     AliError("Null pointer for RecPoints branch");
 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     return -1;
 335   }
 336   AliDebug(1,Form("RecPoints on Layer 1,2 = %d, %d\n",nrpL1,nrpL2));
 337
 338   Double_t a[3]={xbeam,ybeam,0.};
 339   Double_t b[3]={xbeam,ybeam,10.};
 340   AliStrLine zeta(a,b,kTRUE);
 341   static Double_t bField=TMath::Abs(AliTracker::GetBz()/10.); //T
 342   SetMeanPPtSelTracks(bField);
 343
 344   Int_t nolines = 0;
 345   // Loop on modules of layer 1
 346   for(Int_t modul1= firstL1; modul1<=lastL1;modul1++){   // Loop on modules of layer 1
 347     if(!fUseModule[modul1]) continue;
 348     UShort_t ladder=int(modul1/4)+1; // ladders are numbered starting from 1
 349     branch->GetEvent(modul1);
 350     Int_t nrecp1 = itsRec->GetEntries();
 351     static TClonesArray prpl1("AliITSRecPoint",nrecp1);
 352     prpl1.SetOwner();
 353     for(Int_t j=0;j<nrecp1;j++){
 354       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
 355       new(prpl1[j])AliITSRecPoint(*recp);
 356     }
 357     fDetTypeRec->ResetRecPoints();
 358     for(Int_t j=0;j<nrecp1;j++){
 359       if(j>kMaxCluPerMod) continue;
 360       UShort_t idClu1=modul1*kMaxCluPerMod+j;
 361       if(fUsedCluster.TestBitNumber(idClu1)) continue;
 362       AliITSRecPoint *recp1 = (AliITSRecPoint*)prpl1.At(j);
 363       recp1->GetGlobalXYZ(gc1f);
 364       for(Int_t ico=0;ico<3;ico++)gc1[ico]=gc1f[ico];
 365
 366       Double_t phi1 = TMath::ATan2(gc1[1]-ybeam,gc1[0]-xbeam);
 367       if(phi1<0)phi1=2*TMath::Pi()+phi1;
 368       for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){
 369         for(Int_t k=0;k<4;k++){
 370           Int_t ladmod=fLadders[ladder-1]+ladl2;
 371           if(ladmod>AliITSgeomTGeo::GetNLadders(2)) ladmod=ladmod-AliITSgeomTGeo::GetNLadders(2);
 372           Int_t modul2=AliITSgeomTGeo::GetModuleIndex(2,ladmod,k+1);
 373           if(!fUseModule[modul2]) continue;
 374           branch->GetEvent(modul2);
 375           Int_t nrecp2 = itsRec->GetEntries();
 376           for(Int_t j2=0;j2<nrecp2;j2++){
 377             if(j2>kMaxCluPerMod) continue;
 378             UShort_t idClu2=modul2*kMaxCluPerMod+j2;
 379             if(fUsedCluster.TestBitNumber(idClu2)) continue;
 380             AliITSRecPoint *recp2 = (AliITSRecPoint*)itsRec->At(j2);
 381             recp2->GetGlobalXYZ(gc2f);
 382             for(Int_t ico=0;ico<3;ico++)gc2[ico]=gc2f[ico];
 383             Double_t phi2 = TMath::ATan2(gc2[1]-ybeam,gc2[0]-xbeam);
 384             if(phi2<0)phi2=2*TMath::Pi()+phi2;
 385             Double_t diff = TMath::Abs(phi2-phi1);
 386             if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
 387             if(optCuts==0 && diff<fDiffPhiMax){
 388               Double_t r1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]);
 389               Double_t zc1=gc1[2];
 390               Double_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
 391               Double_t zc2=gc2[2];
 392               Double_t zr0=(r2*zc1-r1*zc2)/(r2-r1); //Z @ null radius
 393               fZHisto->Fill(zr0);
 394             }
 395             if(diff>deltaPhi)continue;
 396             AliStrLine line(gc1,gc2,kTRUE);
 397             Double_t cp[3];
 398             Int_t retcode = line.Cross(&zeta,cp);
 399             if(retcode<0)continue;
 400             Double_t dca = line.GetDCA(&zeta);
 401             if(dca<0.) continue;
 402             if(dca>deltaR)continue;
 403             Double_t deltaZ=cp[2]-zvert;
 404             if(TMath::Abs(deltaZ)>dZmax)continue;
 405
 406             if(nolines == 0){
 407               if(fLines.GetEntriesFast()>0)fLines.Clear("C");
 408             }
 409             Float_t cov[6];
 410             recp2->GetGlobalCov(cov);
 411
 412
 413             Double_t rad1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]);
 414             Double_t rad2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
 415             Double_t factor=(rad1+rad2)/(rad2-rad1); //factor to account for error on tracklet direction
 416
 417             Double_t curvErr=0;
 418             if(bField>0.00001){
 419               Double_t curvRadius=fMeanPtSelTrk/(0.3*bField)*100; //cm
 420               Double_t dRad=TMath::Sqrt((gc1[0]-gc2[0])*(gc1[0]-gc2[0])+(gc1[1]-gc2[1])*(gc1[1]-gc2[1]));
 421               Double_t aux=dRad/2.+rad1;
 422               curvErr=TMath::Sqrt(curvRadius*curvRadius-dRad*dRad/4.)-TMath::Sqrt(curvRadius*curvRadius-aux*aux); //cm
 423             }
 424             Double_t sigmasq[3];
 425             sigmasq[0]=(cov[0]+curvErr*curvErr/2.)*factor*factor;
 426             sigmasq[1]=(cov[3]+curvErr*curvErr/2.)*factor*factor;
 427             sigmasq[2]=cov[5]*factor*factor;
 428
 429             // Multiple scattering
 430             Double_t pOverMass=fMeanPSelTrk/0.140;
 431             Double_t beta2=pOverMass*pOverMass/(1+pOverMass*pOverMass);
 432             Double_t p2=fMeanPSelTrk*fMeanPSelTrk;
 433             Double_t rBP=GetPipeRadius();
 434             Double_t dBP=0.08/35.3; // 800 um of Be
 435             Double_t dL1=0.01; //approx. 1% of radiation length
 436             Double_t theta2BP=14.1*14.1/(beta2*p2*1e6)*dBP;
 437             Double_t theta2L1=14.1*14.1/(beta2*p2*1e6)*dL1;
 438             Double_t rtantheta1=(rad2-rad1)*TMath::Tan(TMath::Sqrt(theta2L1));
 439             Double_t rtanthetaBP=(rad1-rBP)*TMath::Tan(TMath::Sqrt(theta2BP));
 440             for(Int_t ico=0; ico<3;ico++){
 441               sigmasq[ico]+=rtantheta1*rtantheta1*factor*factor/3.;
 442               sigmasq[ico]+=rtanthetaBP*rtanthetaBP*factor*factor/3.;
 443             }
 444             Double_t wmat[9]={1.,0.,0.,0.,1.,0.,0.,0.,1.};
 445             if(sigmasq[0]!=0.) wmat[0]=1./sigmasq[0];
 446             if(sigmasq[1]!=0.) wmat[4]=1./sigmasq[1];
 447             if(sigmasq[2]!=0.) wmat[8]=1./sigmasq[2];
 448             new(fLines[nolines++])AliStrLine(gc1,sigmasq,wmat,gc2,kTRUE,idClu1,idClu2);
 449
 450           }
 451           fDetTypeRec->ResetRecPoints();
 452         }
 453       }
 454     }
 455     prpl1.Clear();
 456   }
 457   if(nolines == 0)return -2;
 458   return nolines;
 459 }
 460
 461 //______________________________________________________________________
 462 Int_t  AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){
 463   // Finds the 3D vertex information using tracklets
 464   Int_t retcode = -1;
 465
 466   Double_t xbeam=GetNominalPos()[0];
 467   Double_t ybeam=GetNominalPos()[1];
 468   Double_t zvert=0.;
 469   Double_t deltaR=fCoarseMaxRCut;
 470   if(fPileupAlgo == 2) {
 471     deltaR=fMaxRCut;
 472     if(optCuts!=0)AliWarning(Form("fPileupAlgo=2 AND optCuts=%d. It should be 0",optCuts));
 473   }
 474   Double_t dZmax=fZCutDiamond;
 475   if(optCuts==1){
 476     xbeam=fVert3D.GetXv();
 477     ybeam=fVert3D.GetYv();
 478     zvert=fVert3D.GetZv();
 479     deltaR=fMaxRCut;
 480     dZmax=fMaxZCut;
 481   }else if(optCuts==2){
 482     xbeam=fVert3D.GetXv();
 483     ybeam=fVert3D.GetYv();
 484     deltaR=fMaxRCut;
 485   }
 486
 487   Double_t rl=-fCoarseMaxRCut;
 488   Double_t rh=fCoarseMaxRCut;
 489   Double_t zl=-fZCutDiamond;
 490   Double_t zh=fZCutDiamond;
 491   Int_t nbr=(Int_t)((rh-rl)/fBinSizeR+0.0001);
 492   Int_t nbz=(Int_t)((zh-zl)/fBinSizeZ+0.0001);
 493   Int_t nbrcs=(Int_t)((rh-rl)/(fBinSizeR*2.)+0.0001);
 494   Int_t nbzcs=(Int_t)((zh-zl)/(fBinSizeZ*2.)+0.0001);
 495
 496   TH3F *h3d = NULL;
 497   TH3F *h3dcs = NULL;
 498   if(fPileupAlgo !=2){
 499     h3d = new TH3F("h3d","xyz distribution",nbr,rl,rh,nbr,rl,rh,nbz,zl,zh);
 500     h3dcs = new TH3F("h3dcs","xyz distribution",nbrcs,rl,rh,nbrcs,rl,rh,nbzcs,zl,zh);
 501   }
 502   // cleanup of the TCLonesArray of tracklets (i.e. fakes are removed)
 503   Int_t *validate = new Int_t [fLines.GetEntriesFast()];
 504   for(Int_t i=0; i<fLines.GetEntriesFast();i++)validate[i]=0;
 505   for(Int_t i=0; i<fLines.GetEntriesFast()-1;i++){
 506     AliStrLine *l1 = (AliStrLine*)fLines.At(i);
 507     for(Int_t j=i+1;j<fLines.GetEntriesFast();j++){
 508       AliStrLine *l2 = (AliStrLine*)fLines.At(j);
 509       Double_t dca=l1->GetDCA(l2);
 510       if(dca > fDCAcut || dca<0.00001) continue;
 511       Double_t point[3];
 512       Int_t retc = l1->Cross(l2,point);
 513       if(retc<0)continue;
 514       Double_t deltaZ=point[2]-zvert;
 515      if(TMath::Abs(deltaZ)>dZmax)continue;
 516       Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
 517       if(rad>fCoarseMaxRCut)continue;
 518       Double_t deltaX=point[0]-xbeam;
 519       Double_t deltaY=point[1]-ybeam;
 520       Double_t raddist=TMath::Sqrt(deltaX*deltaX+deltaY*deltaY);
 521       if(raddist>deltaR)continue;
 522       validate[i]=1;
 523       validate[j]=1;
 524       if(fPileupAlgo != 2){
 525         h3d->Fill(point[0],point[1],point[2]);
 526         h3dcs->Fill(point[0],point[1],point[2]);
 527       }
 528     }
 529   }
 530
 531
 532
 533   Int_t numbtracklets=0;
 534   for(Int_t i=0; i<fLines.GetEntriesFast();i++)if(validate[i]>=1)numbtracklets++;
 535   if(numbtracklets<2){
 536     delete [] validate;
 537     if(fPileupAlgo != 2){
 538       delete h3d;
 539       delete h3dcs;
 540     }
 541     return retcode;
 542   }
 543
 544   for(Int_t i=0; i<fLines.GetEntriesFast();i++){
 545     if(validate[i]<1)fLines.RemoveAt(i);
 546   }
 547   fLines.Compress();
 548   AliDebug(1,Form("Number of tracklets (after compress)%d ",fLines.GetEntriesFast()));
 549   delete [] validate;
 550
 551   // Exit here if Pileup Algorithm 2 has been chosen
 552   if(fPileupAlgo == 2)return 0;
 553   //
 554
 555
 556   //        Find peaks in histos
 557
 558   Double_t peak[3]={0.,0.,0.};
 559   Int_t ntrkl,ntimes;
 560   FindPeaks(h3d,peak,ntrkl,ntimes);
 561   delete h3d;
 562   Double_t binsizer=(rh-rl)/nbr;
 563   Double_t binsizez=(zh-zl)/nbz;
 564   if(optCuts==0 && (ntrkl<=2 || ntimes>1)){
 565     ntrkl=0;
 566     ntimes=0;
 567     FindPeaks(h3dcs,peak,ntrkl,ntimes);
 568     binsizer=(rh-rl)/nbrcs;
 569     binsizez=(zh-zl)/nbzcs;
 570     if(ntrkl==1 || ntimes>1){delete h3dcs; return retcode;}
 571   }
 572   delete h3dcs;
 573
 574   //         Second selection loop
 575
 576   Double_t bs=(binsizer+binsizez)/2.;
 577   for(Int_t i=0; i<fLines.GetEntriesFast();i++){
 578     AliStrLine *l1 = (AliStrLine*)fLines.At(i);
 579     if(l1->GetDistFromPoint(peak)>2.5*bs)fLines.RemoveAt(i);
 580   }
 581   fLines.Compress();
 582   AliDebug(1,Form("Number of tracklets (after 2nd compression) %d",fLines.GetEntriesFast()));
 583
 584   if(fLines.GetEntriesFast()>1){
 585     retcode=0;
 586     //  find a first candidate for the primary vertex
 587     fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
 588     // make a further selection on tracklets based on this first candidate
 589     fVert3D.GetXYZ(peak);
 590     AliDebug(1,Form("FIRST V candidate: %f ; %f ; %f",peak[0],peak[1],peak[2]));
 591     for(Int_t i=0; i<fLines.GetEntriesFast();i++){
 592       AliStrLine *l1 = (AliStrLine*)fLines.At(i);
 593       if(l1->GetDistFromPoint(peak)> fDCAcut)fLines.RemoveAt(i);
 594     }
 595     fLines.Compress();
 596     AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines.GetEntriesFast()));
 597     if(fLines.GetEntriesFast()>1){// this new tracklet selection is used
 598       fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
 599     }
 600   }
 601   return retcode;
 602 }
 603
 604 //________________________________________________________
 605 void AliITSVertexer3D::SetMeanPPtSelTracks(Double_t fieldTesla){
 606   // Sets mean values of Pt based on the field
 607   // for P (used in multiple scattering) the most probable value is used
 608   if(TMath::Abs(fieldTesla-0.5)<0.01){
 609     SetMeanPSelTracks(0.375);
 610     SetMeanPtSelTracks(0.630);
 611   }else if(TMath::Abs(fieldTesla-0.4)<0.01){
 612     SetMeanPSelTracks(0.375);
 613     SetMeanPtSelTracks(0.580);
 614   }else if(TMath::Abs(fieldTesla-0.2)<0.01){
 615     SetMeanPSelTracks(0.375);
 616     SetMeanPtSelTracks(0.530);
 617   }else if(fieldTesla<0.00001){
 618     SetMeanPSelTracks(0.375);
 619     SetMeanPtSelTracks(0.230);
 620   }else{
 621     SetMeanPSelTracks();
 622     SetMeanPtSelTracks();
 623   }
 624 }
 625
 626 //________________________________________________________
 627 void AliITSVertexer3D::FindPeaks(TH3F* histo, Double_t *peak, Int_t &nOfTracklets, Int_t &nOfTimes){
 628   // Finds bin with max contents in 3D histo of tracket intersections
 629   TAxis *xax = histo->GetXaxis();
 630   TAxis *yax = histo->GetYaxis();
 631   TAxis *zax = histo->GetZaxis();
 632   peak[0]=0.;
 633   peak[1]=0.;
 634   peak[2]=0.;
 635   nOfTracklets = 0;
 636   nOfTimes=0;
 637   Int_t peakbin[3]={0,0,0};
 638   Int_t peak2bin[3]={-1,-1,-1};
 639   Int_t bc2=-1;
 640   for(Int_t i=xax->GetFirst();i<=xax->GetLast();i++){
 641     Double_t xval = xax->GetBinCenter(i);
 642     for(Int_t j=yax->GetFirst();j<=yax->GetLast();j++){
 643       Double_t yval = yax->GetBinCenter(j);
 644       for(Int_t k=zax->GetFirst();k<=zax->GetLast();k++){
 645         Double_t zval = zax->GetBinCenter(k);
 646         Int_t bc =(Int_t)histo->GetBinContent(i,j,k);
 647         if(bc==0) continue;
 648         if(bc>nOfTracklets){
 649           nOfTracklets=bc;
 650           peak[2] = zval;
 651           peak[1] = yval;
 652           peak[0] = xval;
 653           peakbin[2] = k;
 654           peakbin[1] = j;
 655           peakbin[0] = i;
 656           peak2bin[2] = -1;
 657           peak2bin[1] = -1;
 658           peak2bin[0] = -1;
 659           bc2=-1;
 660           nOfTimes = 1;
 661         }else if(bc==nOfTracklets){
 662           if(TMath::Abs(i-peakbin[0])<=1 && TMath::Abs(j-peakbin[1])<=1 && TMath::Abs(k-peakbin[2])<=1){
 663             peak2bin[2] = k;
 664             peak2bin[1] = j;
 665             peak2bin[0] = i;
 666             bc2=bc;
 667             nOfTimes = 1;
 668           }else{
 669             nOfTimes++;
 670           }
 671         }
 672       }
 673     }
 674   }
 675   if(peak2bin[0]>=-1 && bc2!=-1){ // two contiguous peak-cells with same contents
 676     peak[0]=0.5*(xax->GetBinCenter(peakbin[0])+xax->GetBinCenter(peak2bin[0]));
 677     peak[1]=0.5*(yax->GetBinCenter(peakbin[1])+yax->GetBinCenter(peak2bin[1]));
 678     peak[2]=0.5*(zax->GetBinCenter(peakbin[2])+zax->GetBinCenter(peak2bin[2]));
 679     nOfTracklets+=bc2;
 680     nOfTimes=1;
 681   }
 682 }
 683 //________________________________________________________
 684 void AliITSVertexer3D::MarkUsedClusters(){
 685   // Mark clusters of tracklets used in vertex claulation
 686   for(Int_t i=0; i<fLines.GetEntriesFast();i++){
 687     AliStrLine *lin = (AliStrLine*)fLines.At(i);
 688     Int_t idClu1=lin->GetIdPoint(0);
 689     Int_t idClu2=lin->GetIdPoint(1);
 690     fUsedCluster.SetBitNumber(idClu1);
 691     fUsedCluster.SetBitNumber(idClu2);
 692   }
 693 }
 694 //________________________________________________________
 695 Int_t AliITSVertexer3D::RemoveTracklets(){
 696   // Remove trackelts close to first found vertex
 697   Double_t vert[3]={fVert3D.GetXv(),fVert3D.GetYv(),fVert3D.GetZv()};
 698   Int_t nRemoved=0;
 699   for(Int_t i=0; i<fLines.GetEntriesFast();i++){
 700     AliStrLine *lin = (AliStrLine*)fLines.At(i);
 701     if(lin->GetDistFromPoint(vert)<fDCAforPileup){
 702       Int_t idClu1=lin->GetIdPoint(0);
 703       Int_t idClu2=lin->GetIdPoint(1);
 704       fUsedCluster.SetBitNumber(idClu1);
 705       fUsedCluster.SetBitNumber(idClu2);
 706       fLines.RemoveAt(i);
 707       ++nRemoved;
 708     }
 709   }
 710   fLines.Compress();
 711   return nRemoved;
 712 }
 713 //________________________________________________________
 714 void AliITSVertexer3D::FindOther3DVertices(TTree *itsClusterTree){
 715   // pileup identification based on 3D vertexing with not used clusters
 716   MarkUsedClusters();
 717   fLines.Clear("C");
 718   Int_t nolines = FindTracklets(itsClusterTree,2);
 719   if(nolines>=2){
 720     Int_t nr=RemoveTracklets();
 721     nolines-=nr;
 722     if(nolines>=2){
 723       Int_t rc=Prepare3DVertex(2);
 724       if(rc==0){
 725         fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
 726         if(fVert3D.GetNContributors()>=fMinTrackletsForPilup){
 727           fIsPileup=kTRUE;
 728           fNoVertices=2;
 729           fVertArray = new AliESDVertex[2];
 730           fVertArray[0]=(*fCurrentVertex);
 731           fVertArray[1]=fVert3D;
 732           fZpuv=fVert3D.GetZv();
 733           fNTrpuv=fVert3D.GetNContributors();
 734         }
 735       }
 736     }
 737   }
 738 }
 739 //______________________________________________________________________
 740 void AliITSVertexer3D::PileupFromZ(){
 741   // Calls the pileup algorithm of ALiITSVertexerZ
 742   Int_t binmin, binmax;
 743   Int_t nPeaks=AliITSVertexerZ::GetPeakRegion(fZHisto,binmin,binmax);
 744   if(nPeaks==2)AliWarning("2 peaks found");
 745   Int_t firstPeakCont=0;
 746   Double_t firstPeakPos=0.;
 747   for(Int_t i=binmin-1;i<=binmax+1;i++){
 748     firstPeakCont+=static_cast<Int_t>(fZHisto->GetBinContent(i));
 749     firstPeakPos+=fZHisto->GetBinContent(i)*fZHisto->GetBinCenter(i);
 750   }
 751   if(firstPeakCont>0){
 752     firstPeakPos/=firstPeakCont;
 753     Int_t ncontr2=0;
 754     if(firstPeakCont>fMinTrackletsForPilup){
 755       Float_t secPeakPos;
 756       ncontr2=AliITSVertexerZ::FindSecondPeak(fZHisto,binmin,binmax,secPeakPos);
 757       if(ncontr2>=fMinTrackletsForPilup){
 758         fIsPileup=kTRUE;
 759         fNoVertices=2;
 760         AliESDVertex secondVert(secPeakPos,0.1,ncontr2);
 761         fVertArray = new AliESDVertex[2];
 762         fVertArray[0]=(*fCurrentVertex);
 763         fVertArray[1]=secondVert;
 764         fZpuv=secPeakPos;
 765         fNTrpuv=ncontr2;
 766       }
 767     }
 768   }
 769 }
 770 //________________________________________________________
 771 void AliITSVertexer3D::PrintStatus() const {
 772   // Print current status
 773   printf("=======================================================\n");
 774   printf("Loose cut on Delta Phi %f\n",fCoarseDiffPhiCut);
 775   printf("Cut on tracklet DCA to Z axis %f\n",fCoarseMaxRCut);
 776   printf("Cut on tracklet DCA to beam axis %f\n",fMaxRCut);
 777   printf("Cut on diamond (Z) %f\n",fZCutDiamond);
 778   printf("Cut on DCA - tracklet to tracklet and to vertex %f\n",fDCAcut);
 779   printf("Max Phi difference: %f\n",fDiffPhiMax);
 780   printf("Pileup algo: %d\n",fPileupAlgo);
 781   printf("Min DCA to 1st vetrtex for pileup: %f\n",fDCAforPileup);
 782   printf("=======================================================\n");
 783 }