[u/mrichter/AliRoot.git] / ITS / AliITSVertexer3D.cxx

/**************************************************************************
 * Copyright(c) 2006-2008, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
#include <TTree.h>
#include "AliESDVertex.h"
#include "AliLog.h"
#include "AliStrLine.h"
#include "AliTracker.h"
#include "AliITSDetTypeRec.h"
#include "AliITSRecPoint.h"
#include "AliITSgeomTGeo.h"
#include "AliVertexerTracks.h"
#include "AliITSVertexer3D.h"
#include "AliITSVertexerZ.h"
/////////////////////////////////////////////////////////////////
// this class implements a method to determine
// the 3 coordinates of the primary vertex
// for p-p collisions 
// It can be used successfully with Pb-Pb collisions
////////////////////////////////////////////////////////////////

ClassImp(AliITSVertexer3D)

/* $Id$ */

//______________________________________________________________________
AliITSVertexer3D::AliITSVertexer3D():AliITSVertexer(),
fLines("AliStrLine",1000),
fVert3D(),
fCoarseDiffPhiCut(0.),
fCoarseMaxRCut(0.),
fMaxRCut(0.),
fZCutDiamond(0.),
fMaxZCut(0.),
fDCAcut(0.),
fDiffPhiMax(0.),
fMeanPSelTrk(0.),
fMeanPtSelTrk(0.)
{
  // Default constructor
  SetCoarseDiffPhiCut();
  SetCoarseMaxRCut();
  SetMaxRCut();
  SetZCutDiamond();
  SetMaxZCut();
  SetDCACut();
  SetDiffPhiMax();
  SetMeanPSelTracks();
  SetMeanPtSelTracks();
}

//______________________________________________________________________
AliITSVertexer3D::~AliITSVertexer3D() {
  // Destructor
  fLines.Clear("C");
}

//______________________________________________________________________
void AliITSVertexer3D::ResetVert3D(){
  //
  fVert3D.SetXv(0.);
  fVert3D.SetYv(0.);
  fVert3D.SetZv(0.);
  fVert3D.SetDispersion(0.);
  fVert3D.SetNContributors(0);
}
//______________________________________________________________________
AliESDVertex* AliITSVertexer3D::FindVertexForCurrentEvent(TTree *itsClusterTree){
  // Defines the AliESDVertex for the current event
  ResetVert3D();
  AliDebug(1,"FindVertexForCurrentEvent - 3D - PROCESSING NEXT EVENT");
  fLines.Clear();

  Int_t nolines = FindTracklets(itsClusterTree,0);
  fCurrentVertex = 0;
  if(nolines>=2){
    Int_t rc=Prepare3DVertex(0);
    if(rc==0) fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
  /*  uncomment to debug
      printf("Vertex found in first iteration:\n");
      fVert3D.Print();
      printf("Start second iteration\n");
  end of debug lines  */
    if(fVert3D.GetNContributors()>0){
      fLines.Clear("C");
      nolines = FindTracklets(itsClusterTree,1);
      if(nolines>=2){
	rc=Prepare3DVertex(1);
	if(rc==0) fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
      }
    }
  /*  uncomment to debug 
      printf("Vertex found in second iteration:\n");
      fVert3D.Print();
      end of debug lines  */ 
 
    Float_t vRadius=TMath::Sqrt(fVert3D.GetXv()*fVert3D.GetXv()+fVert3D.GetYv()*fVert3D.GetYv());
    if(vRadius<GetPipeRadius() && fVert3D.GetNContributors()>0){
      Double_t position[3]={fVert3D.GetXv(),fVert3D.GetYv(),fVert3D.GetZv()};
      Double_t covmatrix[6];
      fVert3D.GetCovMatrix(covmatrix);
      Double_t chi2=99999.;
      Int_t    nContr=fVert3D.GetNContributors();
      fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr);    
      fCurrentVertex->SetTitle("vertexer: 3D");
      fCurrentVertex->SetName("SPDVertex3D");
      fCurrentVertex->SetDispersion(fVert3D.GetDispersion());
    }
  }
  if(!fCurrentVertex){
    AliITSVertexerZ vertz(GetNominalPos()[0],GetNominalPos()[1]);
    vertz.SetDetTypeRec(GetDetTypeRec());
    AliDebug(1,"Call Vertexer Z\n");
    vertz.SetLowLimit(-fZCutDiamond);
    vertz.SetHighLimit(fZCutDiamond);
    AliESDVertex* vtxz = vertz.FindVertexForCurrentEvent(itsClusterTree);
    if(vtxz){
      Double_t position[3]={GetNominalPos()[0],GetNominalPos()[1],vtxz->GetZv()};
      Double_t covmatrix[6];
      vtxz->GetCovMatrix(covmatrix);
      Double_t chi2=99999.;
      Int_t    nContr=vtxz->GetNContributors();
      fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr);    
      fCurrentVertex->SetTitle("vertexer: Z");
      fCurrentVertex->SetName("SPDVertexZ");
      delete vtxz;
      //      printf("Vertexer Z success\n");
    }

  }
  FindMultiplicity(itsClusterTree);
  return fCurrentVertex;
}  

//______________________________________________________________________
Int_t AliITSVertexer3D::FindTracklets(TTree *itsClusterTree, Int_t optCuts){
  // All the possible combinations between recpoints on layer 1and 2 are
  // considered. Straight lines (=tracklets)are formed. 
  // The tracklets are processed in Prepare3DVertex
  AliITSDetTypeRec detTypeRec;

  TTree *tR = itsClusterTree;
  detTypeRec.SetTreeAddressR(tR);
  TClonesArray *itsRec  = 0;
  // lc1 and gc1 are local and global coordinates for layer 1
  //  Float_t lc1[3]={0.,0.,0.};
  Float_t gc1[3]={0.,0.,0.};
  // lc2 and gc2 are local and global coordinates for layer 2
  //  Float_t lc2[3]={0.,0.,0.};
  Float_t gc2[3]={0.,0.,0.};

  itsRec = detTypeRec.RecPoints();
  TBranch *branch;
  branch = tR->GetBranch("ITSRecPoints");

  // Set values for cuts
  Float_t xbeam=GetNominalPos()[0]; 
  Float_t ybeam=GetNominalPos()[1];
  Float_t zvert=0.;
  Float_t deltaPhi=fCoarseDiffPhiCut;
  Float_t deltaR=fCoarseMaxRCut;
  Float_t dZmax=fZCutDiamond;
  if(optCuts){
    xbeam=fVert3D.GetXv();
    ybeam=fVert3D.GetYv();
    zvert=fVert3D.GetZv();
    deltaPhi = fDiffPhiMax; 
    deltaR=fMaxRCut;
    dZmax=fMaxZCut;
  }
  Int_t nrpL1 = 0;    // number of rec points on layer 1
  Int_t nrpL2 = 0;    // number of rec points on layer 2

  // By default irstL1=0 and lastL1=79
  Int_t firstL1 = AliITSgeomTGeo::GetModuleIndex(1,1,1);
  Int_t lastL1 = AliITSgeomTGeo::GetModuleIndex(2,1,1)-1;
  for(Int_t module= firstL1; module<=lastL1;module++){  // count number of recopints on layer 1
    branch->GetEvent(module);
    nrpL1+= itsRec->GetEntries();
    detTypeRec.ResetRecPoints();
  }
  //By default firstL2=80 and lastL2=239
  Int_t firstL2 = AliITSgeomTGeo::GetModuleIndex(2,1,1);
  Int_t lastL2 = AliITSgeomTGeo::GetModuleIndex(3,1,1)-1;
  for(Int_t module= firstL2; module<=lastL2;module++){  // count number of recopints on layer 2
    branch->GetEvent(module);
    nrpL2+= itsRec->GetEntries();
    detTypeRec.ResetRecPoints();
  }
  if(nrpL1 == 0 || nrpL2 == 0){
    return -1;
  }
  AliDebug(1,Form("RecPoints on Layer 1,2 = %d, %d\n",nrpL1,nrpL2));

  Double_t a[3]={xbeam,ybeam,0.}; 
  Double_t b[3]={xbeam,ybeam,10.};
  AliStrLine zeta(a,b,kTRUE);
  Float_t bField=AliTracker::GetBz()/10.; //T
  SetMeanPPtSelTracks(bField);

  Int_t nolines = 0;
  // Loop on modules of layer 1
  for(Int_t modul1= firstL1; modul1<=lastL1;modul1++){   // Loop on modules of layer 1
    if(!fUseModule[modul1]) continue;
    UShort_t ladder=int(modul1/4)+1; // ladders are numbered starting from 1
    branch->GetEvent(modul1);
    Int_t nrecp1 = itsRec->GetEntries();
    static TClonesArray prpl1("AliITSRecPoint",nrecp1);
    prpl1.SetOwner();
    for(Int_t j=0;j<nrecp1;j++){
      AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
      new(prpl1[j])AliITSRecPoint(*recp);
    }
    detTypeRec.ResetRecPoints();
    for(Int_t j=0;j<nrecp1;j++){
      AliITSRecPoint *recp1 = (AliITSRecPoint*)prpl1.At(j);
      // Local coordinates of this recpoint
      /*
      lc[0]=recp1->GetDetLocalX();
      lc[2]=recp1->GetDetLocalZ();
      */
      recp1->GetGlobalXYZ(gc1);
      Double_t phi1 = TMath::ATan2(gc1[1]-ybeam,gc1[0]-xbeam);
      if(phi1<0)phi1=2*TMath::Pi()+phi1;
      for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){
	for(Int_t k=0;k<4;k++){
	  Int_t ladmod=fLadders[ladder-1]+ladl2;
 	  if(ladmod>AliITSgeomTGeo::GetNLadders(2)) ladmod=ladmod-AliITSgeomTGeo::GetNLadders(2);
	  Int_t modul2=AliITSgeomTGeo::GetModuleIndex(2,ladmod,k+1);
	  if(!fUseModule[modul2]) continue;
	  branch->GetEvent(modul2);
	  Int_t nrecp2 = itsRec->GetEntries();
	  for(Int_t j2=0;j2<nrecp2;j2++){
	    AliITSRecPoint *recp2 = (AliITSRecPoint*)itsRec->At(j2);
	    /*
	    lc2[0]=recp2->GetDetLocalX();
	    lc2[2]=recp2->GetDetLocalZ();
	    */
	    recp2->GetGlobalXYZ(gc2);
	    Double_t phi2 = TMath::ATan2(gc2[1]-ybeam,gc2[0]-xbeam);
	    if(phi2<0)phi2=2*TMath::Pi()+phi2;
	    Double_t diff = TMath::Abs(phi2-phi1); 
	    if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff; 
	    if(diff>deltaPhi)continue;
	    AliStrLine line(gc1,gc2,kTRUE);
	    Double_t cp[3];
	    Int_t retcode = line.Cross(&zeta,cp);
	    if(retcode<0)continue;
	    Double_t dca = line.GetDCA(&zeta);
	    if(dca<0.) continue;
	    if(dca>deltaR)continue;
	    Double_t deltaZ=cp[2]-zvert;
	    if(TMath::Abs(deltaZ)>dZmax)continue;

	    if(nolines == 0){
	      if(fLines.GetEntriesFast()>0)fLines.Clear("C");
	    }
	    Float_t cov[6];
	    recp2->GetGlobalCov(cov);

	    
	    Float_t rad1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]);
	    Float_t rad2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
	    Float_t factor=(rad1+rad2)/(rad2-rad1); //factor to account for error on tracklet direction 

	    Float_t curvErr=0;
	    if(bField>0.00001){
	      Float_t curvRadius=fMeanPtSelTrk/(0.3*bField)*100; //cm 
	      Float_t dRad=TMath::Sqrt(TMath::Power((gc1[0]-gc2[0]),2)+TMath::Power((gc1[1]-gc2[1]),2));
	      Float_t aux=dRad/2.+rad1;
	      curvErr=TMath::Sqrt(curvRadius*curvRadius-dRad*dRad/4.)-TMath::Sqrt(curvRadius*curvRadius-aux*aux); //cm
	    }

	    Float_t sigmasq[3];
	    sigmasq[0]=(cov[0]+curvErr*curvErr/2.)*factor*factor;
	    sigmasq[1]=(cov[3]+curvErr*curvErr/2.)*factor*factor;
	    sigmasq[2]=cov[5]*factor*factor;

	    // Multiple scattering
 	    Float_t beta=1.;
 	    Float_t beta2=beta*beta;
 	    Float_t p2=fMeanPSelTrk*fMeanPSelTrk;
 	    Float_t rBP=GetPipeRadius();
 	    Float_t dBP=0.08/35.3; // 800 um of Be
 	    Float_t dL1=0.01; //approx. 1% of radiation length  
 	    Float_t theta2BP=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(dBP);
 	    Float_t theta2L1=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(dL1);
 	    Float_t thetaBP=TMath::Sqrt(theta2BP);
 	    Float_t thetaL1=TMath::Sqrt(theta2L1);
//  	    Float_t geomfac[3];
//  	    geomfac[0]=sin(phi1)*sin(phi1);
//  	    geomfac[1]=cos(phi1)*cos(phi1);
//  	    Float_t tgth=(gc2[2]-gc1[2])/(rad2-rad1);	    
//  	    geomfac[2]=1+tgth*tgth;
 	    for(Int_t ico=0; ico<3;ico++){    
// 	      printf("Error on coord. %d due to cov matrix+curvErr=%f\n",ico,sigmasq[ico]);
// 	      //	      sigmasq[ico]+=rad1*rad1*geomfac[ico]*theta2L1/2; // multiple scattering in layer 1
// 	      //  sigmasq[ico]+=rBP*rBP*geomfac[ico]*theta2BP/2; // multiple scattering in beam pipe
 	      sigmasq[ico]+=TMath::Power(rad1*TMath::Tan(thetaL1),2)/3.;
 	      sigmasq[ico]+=TMath::Power(rBP*TMath::Tan(thetaBP),2)/3.;

// 	      printf("Multipl. scatt. contr %d = %f (LAY1), %f (BP)\n",ico,rad1*rad1*geomfac[ico]*theta2L1/2,rBP*rBP*geomfac[ico]*theta2BP/2);
// 	      printf("Total error on coord %d = %f\n",ico,sigmasq[ico]);
 	    }
	    Float_t wmat[9]={1.,0.,0.,0.,1.,0.,0.,0.,1.};
	    if(sigmasq[0]!=0.) wmat[0]=1./sigmasq[0];
	    if(sigmasq[1]!=0.) wmat[4]=1./sigmasq[1];
	    if(sigmasq[2]!=0.) wmat[8]=1./sigmasq[2];
	    new(fLines[nolines++])AliStrLine(gc1,sigmasq,wmat,gc2,kTRUE);

	  }
	  detTypeRec.ResetRecPoints();
	}
      }
    }
    prpl1.Clear();
  }
  if(nolines == 0)return -2;
  return nolines;
}

//______________________________________________________________________
Int_t  AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){
  // Finds the 3D vertex information using tracklets
  Int_t retcode = -1;

  Float_t xbeam=GetNominalPos()[0];
  Float_t ybeam=GetNominalPos()[1];
  Float_t zvert=0.;
  Float_t deltaR=fCoarseMaxRCut;
  Float_t dZmax=fZCutDiamond;
  if(optCuts){
    xbeam=fVert3D.GetXv();
    ybeam=fVert3D.GetYv();
    zvert=fVert3D.GetZv();
    deltaR=fMaxRCut;
    dZmax=fMaxZCut;
  }

  Int_t nbr=50;
  Float_t rl=-fCoarseMaxRCut;
  Float_t rh=fCoarseMaxRCut;
  Int_t nbz=100;
  Float_t zl=-fZCutDiamond;
  Float_t zh=fZCutDiamond;
  Float_t binsizer=(rh-rl)/nbr;
  Float_t binsizez=(zh-zl)/nbz;
  TH3F *h3d = new TH3F("h3d","xyz distribution",nbr,rl,rh,nbr,rl,rh,nbz,zl,zh);
  Int_t nbrcs=25;
  Int_t nbzcs=50;
  TH3F *h3dcs = new TH3F("h3dcs","xyz distribution",nbrcs,rl,rh,nbrcs,rl,rh,nbzcs,zl,zh);

  // cleanup of the TCLonesArray of tracklets (i.e. fakes are removed)
  Int_t *validate = new Int_t [fLines.GetEntriesFast()];
  for(Int_t i=0; i<fLines.GetEntriesFast();i++)validate[i]=0;
  for(Int_t i=0; i<fLines.GetEntriesFast()-1;i++){
    AliStrLine *l1 = (AliStrLine*)fLines.At(i);
    for(Int_t j=i+1;j<fLines.GetEntriesFast();j++){
      AliStrLine *l2 = (AliStrLine*)fLines.At(j);
      Float_t dca=l1->GetDCA(l2);
      if(dca > fDCAcut || dca<0.00001) continue;
      Double_t point[3];
      Int_t retc = l1->Cross(l2,point);
      if(retc<0)continue;
      Double_t deltaZ=point[2]-zvert;
     if(TMath::Abs(deltaZ)>dZmax)continue;
      Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
      if(rad>fCoarseMaxRCut)continue;
      Double_t deltaX=point[0]-xbeam;
      Double_t deltaY=point[1]-ybeam;
      Double_t raddist=TMath::Sqrt(deltaX*deltaX+deltaY*deltaY);
      if(raddist>deltaR)continue;
      validate[i]=1;
      validate[j]=1;
      h3d->Fill(point[0],point[1],point[2]);
      h3dcs->Fill(point[0],point[1],point[2]);
    }
  }


  Int_t numbtracklets=0;
  for(Int_t i=0; i<fLines.GetEntriesFast();i++)if(validate[i]>=1)numbtracklets++;
  if(numbtracklets<2){delete [] validate; delete h3d; delete h3dcs; return retcode; }

  for(Int_t i=0; i<fLines.GetEntriesFast();i++){
    if(validate[i]<1)fLines.RemoveAt(i);
  }
  fLines.Compress();
  AliDebug(1,Form("Number of tracklets (after compress)%d ",fLines.GetEntriesFast()));
  delete [] validate;

  //        Find peaks in histos

  Double_t peak[3]={0.,0.,0.};
  Int_t ntrkl,ntimes;
  FindPeaks(h3d,peak,ntrkl,ntimes);  
  delete h3d;

  if(optCuts==0 && ntrkl<=2){
    ntrkl=0;
    ntimes=0;
    FindPeaks(h3dcs,peak,ntrkl,ntimes);  
    binsizer=(rh-rl)/nbrcs;
    binsizez=(zh-zl)/nbzcs;
    if(ntrkl==1 || ntimes>1){delete h3dcs; return retcode;}
  }
  delete h3dcs;

  //         Second selection loop

  Float_t bs=(binsizer+binsizez)/2.;
  for(Int_t i=0; i<fLines.GetEntriesFast();i++){
    AliStrLine *l1 = (AliStrLine*)fLines.At(i);
    if(l1->GetDistFromPoint(peak)>2.5*bs)fLines.RemoveAt(i);
  }
  fLines.Compress();
  AliDebug(1,Form("Number of tracklets (after 2nd compression) %d",fLines.GetEntriesFast()));

  if(fLines.GetEntriesFast()>1){
    //  find a first candidate for the primary vertex
    fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); 
    // make a further selection on tracklets based on this first candidate
    fVert3D.GetXYZ(peak);
    AliDebug(1,Form("FIRST V candidate: %f ; %f ; %f",peak[0],peak[1],peak[2]));
    for(Int_t i=0; i<fLines.GetEntriesFast();i++){
      AliStrLine *l1 = (AliStrLine*)fLines.At(i);
      if(l1->GetDistFromPoint(peak)> fDCAcut)fLines.RemoveAt(i);
    }
    fLines.Compress();
    AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines.GetEntriesFast()));
    if(fLines.GetEntriesFast()>1) retcode=0; // this new tracklet selection is used
    else retcode =1; // the previous tracklet selection will be used
  }
  else {
    retcode = 0;
  }
  return retcode;  
}

//________________________________________________________
void AliITSVertexer3D::SetMeanPPtSelTracks(Float_t fieldTesla){
  // Sets mean values of P and Pt based on the field
  if(TMath::Abs(fieldTesla-0.5)<0.01){
    SetMeanPSelTracks(0.885);
    SetMeanPtSelTracks(0.630);
  }else if(TMath::Abs(fieldTesla-0.4)<0.01){
    SetMeanPSelTracks(0.805);
    SetMeanPtSelTracks(0.580);
  }else if(TMath::Abs(fieldTesla-0.2)<0.01){
    SetMeanPSelTracks(0.740);
    SetMeanPtSelTracks(0.530);
  }else if(fieldTesla<0.00001){
    SetMeanPSelTracks(0.730);
    SetMeanPtSelTracks(0.510);
  }else{
    SetMeanPSelTracks();
    SetMeanPtSelTracks();
  }
}

//________________________________________________________
void AliITSVertexer3D::FindPeaks(TH3F* histo, Double_t *peak, Int_t &nOfTracklets, Int_t &nOfTimes){
  // Finds bin with max contents in 3D histo of tracket intersections
  TAxis *xax = histo->GetXaxis();  
  TAxis *yax = histo->GetYaxis();
  TAxis *zax = histo->GetZaxis();
  peak[0]=0.;
  peak[1]=0.;
  peak[2]=0.;
  nOfTracklets = 0;
  nOfTimes=0;
  for(Int_t i=xax->GetFirst();i<=xax->GetLast();i++){
    Float_t xval = xax->GetBinCenter(i);
    for(Int_t j=yax->GetFirst();j<=yax->GetLast();j++){
      Float_t yval = yax->GetBinCenter(j);
      for(Int_t k=zax->GetFirst();k<=zax->GetLast();k++){
	Float_t zval = zax->GetBinCenter(k);
	Int_t bc =(Int_t)histo->GetBinContent(i,j,k);
	if(bc>nOfTracklets){
	  nOfTracklets = bc;
	  peak[2] = zval;
	  peak[1] = yval;
	  peak[0] = xval;
	  nOfTimes = 1;
	}
	if(bc==nOfTracklets){
	  nOfTimes++;
	}
      }
    }
  }
  
}

//________________________________________________________
void AliITSVertexer3D::PrintStatus() const {
  // Print current status
  printf("=======================================================\n");
  printf("Loose cut on Delta Phi %f\n",fCoarseDiffPhiCut);
  printf("Cut on tracklet DCA to Z axis %f\n",fCoarseMaxRCut);
  printf("Cut on tracklet DCA to beam axis %f\n",fMaxRCut);
  printf("Cut on diamond (Z) %f\n",fZCutDiamond);
  printf("Cut on DCA - tracklet to tracklet and to vertex %f\n",fDCAcut);
  printf(" Max Phi difference: %f\n",fDiffPhiMax);
  printf("=======================================================\n");
}
Commit	Line	Data
70c95f95	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	**************************************************************************/
70c95f95	15	#include <TTree.h>
6a0d56b8	16	#include "AliESDVertex.h"
cd706e57	17	#include "AliLog.h"
6a0d56b8	18	#include "AliStrLine.h"
6a0d56b8	19	#include "AliTracker.h"
70c95f95	20	#include "AliITSDetTypeRec.h"
70c95f95	21	#include "AliITSRecPoint.h"
32e63e47	22	#include "AliITSgeomTGeo.h"
6a0d56b8	23	#include "AliVertexerTracks.h"
6a0d56b8	24	#include "AliITSVertexer3D.h"
0b7bac55	25	#include "AliITSVertexerZ.h"
70c95f95	26	/////////////////////////////////////////////////////////////////
	27	// this class implements a method to determine
	28	// the 3 coordinates of the primary vertex
	29	// for p-p collisions
	30	// It can be used successfully with Pb-Pb collisions
	31	////////////////////////////////////////////////////////////////
	32
	33	ClassImp(AliITSVertexer3D)
	34
6a0d56b8	35	/* $Id$ */
6a0d56b8	36
70c95f95	37	//______________________________________________________________________
70c95f95	38	AliITSVertexer3D::AliITSVertexer3D():AliITSVertexer(),
308c2f7c	39	fLines("AliStrLine",1000),
70c95f95	40	fVert3D(),
70c95f95	41	fCoarseDiffPhiCut(0.),
05d1294c	42	fCoarseMaxRCut(0.),
70c95f95	43	fMaxRCut(0.),
70c95f95	44	fZCutDiamond(0.),
05d1294c	45	fMaxZCut(0.),
70c95f95	46	fDCAcut(0.),
6a0d56b8	47	fDiffPhiMax(0.),
	48	fMeanPSelTrk(0.),
	49	fMeanPtSelTrk(0.)
	50	{
70c95f95	51	// Default constructor
70c95f95	52	SetCoarseDiffPhiCut();
05d1294c	53	SetCoarseMaxRCut();
70c95f95	54	SetMaxRCut();
70c95f95	55	SetZCutDiamond();
05d1294c	56	SetMaxZCut();
7203e11a	57	SetDCACut();
70c95f95	58	SetDiffPhiMax();
6a0d56b8	59	SetMeanPSelTracks();
6a0d56b8	60	SetMeanPtSelTracks();
70c95f95	61	}
70c95f95	62
70c95f95	63	//______________________________________________________________________
	64	AliITSVertexer3D::~AliITSVertexer3D() {
	65	// Destructor
f5f6da22	66	fLines.Clear("C");
70c95f95	67	}
70c95f95	68
27167524	69	//______________________________________________________________________
	70	void AliITSVertexer3D::ResetVert3D(){
	71	//
	72	fVert3D.SetXv(0.);
	73	fVert3D.SetYv(0.);
	74	fVert3D.SetZv(0.);
	75	fVert3D.SetDispersion(0.);
	76	fVert3D.SetNContributors(0);
	77	}
70c95f95	78	//______________________________________________________________________
308c2f7c	79	AliESDVertex* AliITSVertexer3D::FindVertexForCurrentEvent(TTree *itsClusterTree){
70c95f95	80	// Defines the AliESDVertex for the current event
27167524	81	ResetVert3D();
308c2f7c	82	AliDebug(1,"FindVertexForCurrentEvent - 3D - PROCESSING NEXT EVENT");
f5f6da22	83	fLines.Clear();
70c95f95	84
308c2f7c	85	Int_t nolines = FindTracklets(itsClusterTree,0);
70c95f95	86	fCurrentVertex = 0;
0b7bac55	87	if(nolines>=2){
	88	Int_t rc=Prepare3DVertex(0);
	89	if(rc==0) fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
27167524	90	/* uncomment to debug
0b7bac55	91	printf("Vertex found in first iteration:\n");
	92	fVert3D.Print();
	93	printf("Start second iteration\n");
27167524	94	end of debug lines */
0b7bac55	95	if(fVert3D.GetNContributors()>0){
	96	fLines.Clear("C");
	97	nolines = FindTracklets(itsClusterTree,1);
	98	if(nolines>=2){
	99	rc=Prepare3DVertex(1);
	100	if(rc==0) fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
	101	}
d35ad08f	102	}
27167524	103	/* uncomment to debug
0b7bac55	104	printf("Vertex found in second iteration:\n");
	105	fVert3D.Print();
	106	end of debug lines */
27167524	107
0b7bac55	108	Float_t vRadius=TMath::Sqrt(fVert3D.GetXv()fVert3D.GetXv()+fVert3D.GetYv()fVert3D.GetYv());
	109	if(vRadius<GetPipeRadius() && fVert3D.GetNContributors()>0){
	110	Double_t position[3]={fVert3D.GetXv(),fVert3D.GetYv(),fVert3D.GetZv()};
	111	Double_t covmatrix[6];
	112	fVert3D.GetCovMatrix(covmatrix);
	113	Double_t chi2=99999.;
	114	Int_t nContr=fVert3D.GetNContributors();
	115	fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr);
	116	fCurrentVertex->SetTitle("vertexer: 3D");
	117	fCurrentVertex->SetName("SPDVertex3D");
	118	fCurrentVertex->SetDispersion(fVert3D.GetDispersion());
	119	}
	120	}
	121	if(!fCurrentVertex){
	122	AliITSVertexerZ vertz(GetNominalPos()[0],GetNominalPos()[1]);
b96ee725	123	vertz.SetDetTypeRec(GetDetTypeRec());
0b7bac55	124	AliDebug(1,"Call Vertexer Z\n");
cc282e54	125	vertz.SetLowLimit(-fZCutDiamond);
cc282e54	126	vertz.SetHighLimit(fZCutDiamond);
0b7bac55	127	AliESDVertex* vtxz = vertz.FindVertexForCurrentEvent(itsClusterTree);
	128	if(vtxz){
	129	Double_t position[3]={GetNominalPos()[0],GetNominalPos()[1],vtxz->GetZv()};
	130	Double_t covmatrix[6];
	131	vtxz->GetCovMatrix(covmatrix);
	132	Double_t chi2=99999.;
	133	Int_t nContr=vtxz->GetNContributors();
	134	fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr);
	135	fCurrentVertex->SetTitle("vertexer: Z");
	136	fCurrentVertex->SetName("SPDVertexZ");
	137	delete vtxz;
	138	// printf("Vertexer Z success\n");
	139	}
	140
70c95f95	141	}
308c2f7c	142	FindMultiplicity(itsClusterTree);
70c95f95	143	return fCurrentVertex;
	144	}
	145
	146	//______________________________________________________________________
308c2f7c	147	Int_t AliITSVertexer3D::FindTracklets(TTree *itsClusterTree, Int_t optCuts){
70c95f95	148	// All the possible combinations between recpoints on layer 1and 2 are
	149	// considered. Straight lines (=tracklets)are formed.
	150	// The tracklets are processed in Prepare3DVertex
70c95f95	151	AliITSDetTypeRec detTypeRec;
70c95f95	152
308c2f7c	153	TTree *tR = itsClusterTree;
70c95f95	154	detTypeRec.SetTreeAddressR(tR);
70c95f95	155	TClonesArray *itsRec = 0;
6a0d56b8	156	// lc1 and gc1 are local and global coordinates for layer 1
	157	// Float_t lc1[3]={0.,0.,0.};
	158	Float_t gc1[3]={0.,0.,0.};
70c95f95	159	// lc2 and gc2 are local and global coordinates for layer 2
6a0d56b8	160	// Float_t lc2[3]={0.,0.,0.};
6a0d56b8	161	Float_t gc2[3]={0.,0.,0.};
70c95f95	162
	163	itsRec = detTypeRec.RecPoints();
	164	TBranch *branch;
	165	branch = tR->GetBranch("ITSRecPoints");
	166
	167	// Set values for cuts
ecce5370	168	Float_t xbeam=GetNominalPos()[0];
ecce5370	169	Float_t ybeam=GetNominalPos()[1];
05d1294c	170	Float_t zvert=0.;
70c95f95	171	Float_t deltaPhi=fCoarseDiffPhiCut;
05d1294c	172	Float_t deltaR=fCoarseMaxRCut;
05d1294c	173	Float_t dZmax=fZCutDiamond;
70c95f95	174	if(optCuts){
27167524	175	xbeam=fVert3D.GetXv();
	176	ybeam=fVert3D.GetYv();
	177	zvert=fVert3D.GetZv();
70c95f95	178	deltaPhi = fDiffPhiMax;
05d1294c	179	deltaR=fMaxRCut;
05d1294c	180	dZmax=fMaxZCut;
70c95f95	181	}
70c95f95	182	Int_t nrpL1 = 0; // number of rec points on layer 1
	183	Int_t nrpL2 = 0; // number of rec points on layer 2
	184
	185	// By default irstL1=0 and lastL1=79
6a0d56b8	186	Int_t firstL1 = AliITSgeomTGeo::GetModuleIndex(1,1,1);
32e63e47	187	Int_t lastL1 = AliITSgeomTGeo::GetModuleIndex(2,1,1)-1;
6a0d56b8	188	for(Int_t module= firstL1; module<=lastL1;module++){ // count number of recopints on layer 1
70c95f95	189	branch->GetEvent(module);
	190	nrpL1+= itsRec->GetEntries();
	191	detTypeRec.ResetRecPoints();
	192	}
6a0d56b8	193	//By default firstL2=80 and lastL2=239
6a0d56b8	194	Int_t firstL2 = AliITSgeomTGeo::GetModuleIndex(2,1,1);
32e63e47	195	Int_t lastL2 = AliITSgeomTGeo::GetModuleIndex(3,1,1)-1;
6a0d56b8	196	for(Int_t module= firstL2; module<=lastL2;module++){ // count number of recopints on layer 2
70c95f95	197	branch->GetEvent(module);
	198	nrpL2+= itsRec->GetEntries();
	199	detTypeRec.ResetRecPoints();
	200	}
	201	if(nrpL1 == 0 \|\| nrpL2 == 0){
70c95f95	202	return -1;
70c95f95	203	}
cd706e57	204	AliDebug(1,Form("RecPoints on Layer 1,2 = %d, %d\n",nrpL1,nrpL2));
70c95f95	205
	206	Double_t a[3]={xbeam,ybeam,0.};
	207	Double_t b[3]={xbeam,ybeam,10.};
	208	AliStrLine zeta(a,b,kTRUE);
6a0d56b8	209	Float_t bField=AliTracker::GetBz()/10.; //T
6a0d56b8	210	SetMeanPPtSelTracks(bField);
70c95f95	211
	212	Int_t nolines = 0;
	213	// Loop on modules of layer 1
6a0d56b8	214	for(Int_t modul1= firstL1; modul1<=lastL1;modul1++){ // Loop on modules of layer 1
8c42830a	215	if(!fUseModule[modul1]) continue;
27167524	216	UShort_t ladder=int(modul1/4)+1; // ladders are numbered starting from 1
70c95f95	217	branch->GetEvent(modul1);
70c95f95	218	Int_t nrecp1 = itsRec->GetEntries();
f5f6da22	219	static TClonesArray prpl1("AliITSRecPoint",nrecp1);
f5f6da22	220	prpl1.SetOwner();
70c95f95	221	for(Int_t j=0;j<nrecp1;j++){
70c95f95	222	AliITSRecPoint recp = (AliITSRecPoint)itsRec->At(j);
f5f6da22	223	new(prpl1[j])AliITSRecPoint(*recp);
70c95f95	224	}
	225	detTypeRec.ResetRecPoints();
	226	for(Int_t j=0;j<nrecp1;j++){
f5f6da22	227	AliITSRecPoint recp1 = (AliITSRecPoint)prpl1.At(j);
70c95f95	228	// Local coordinates of this recpoint
32e63e47	229	/*
6a0d56b8	230	lc[0]=recp1->GetDetLocalX();
6a0d56b8	231	lc[2]=recp1->GetDetLocalZ();
32e63e47	232	*/
6a0d56b8	233	recp1->GetGlobalXYZ(gc1);
6a0d56b8	234	Double_t phi1 = TMath::ATan2(gc1[1]-ybeam,gc1[0]-xbeam);
70c95f95	235	if(phi1<0)phi1=2*TMath::Pi()+phi1;
27167524	236	for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){
	237	for(Int_t k=0;k<4;k++){
	238	Int_t ladmod=fLadders[ladder-1]+ladl2;
32e63e47	239	if(ladmod>AliITSgeomTGeo::GetNLadders(2)) ladmod=ladmod-AliITSgeomTGeo::GetNLadders(2);
32e63e47	240	Int_t modul2=AliITSgeomTGeo::GetModuleIndex(2,ladmod,k+1);
8c42830a	241	if(!fUseModule[modul2]) continue;
27167524	242	branch->GetEvent(modul2);
	243	Int_t nrecp2 = itsRec->GetEntries();
	244	for(Int_t j2=0;j2<nrecp2;j2++){
6a0d56b8	245	AliITSRecPoint recp2 = (AliITSRecPoint)itsRec->At(j2);
32e63e47	246	/*
6a0d56b8	247	lc2[0]=recp2->GetDetLocalX();
6a0d56b8	248	lc2[2]=recp2->GetDetLocalZ();
32e63e47	249	*/
6a0d56b8	250	recp2->GetGlobalXYZ(gc2);
27167524	251	Double_t phi2 = TMath::ATan2(gc2[1]-ybeam,gc2[0]-xbeam);
	252	if(phi2<0)phi2=2*TMath::Pi()+phi2;
	253	Double_t diff = TMath::Abs(phi2-phi1);
	254	if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
	255	if(diff>deltaPhi)continue;
6a0d56b8	256	AliStrLine line(gc1,gc2,kTRUE);
27167524	257	Double_t cp[3];
	258	Int_t retcode = line.Cross(&zeta,cp);
	259	if(retcode<0)continue;
	260	Double_t dca = line.GetDCA(&zeta);
	261	if(dca<0.) continue;
	262	if(dca>deltaR)continue;
	263	Double_t deltaZ=cp[2]-zvert;
	264	if(TMath::Abs(deltaZ)>dZmax)continue;
6a0d56b8	265
6a0d56b8	266	if(nolines == 0){
f5f6da22	267	if(fLines.GetEntriesFast()>0)fLines.Clear("C");
6a0d56b8	268	}
	269	Float_t cov[6];
	270	recp2->GetGlobalCov(cov);
	271
	272
	273	Float_t rad1=TMath::Sqrt(gc1[0]gc1[0]+gc1[1]gc1[1]);
	274	Float_t rad2=TMath::Sqrt(gc2[0]gc2[0]+gc2[1]gc2[1]);
	275	Float_t factor=(rad1+rad2)/(rad2-rad1); //factor to account for error on tracklet direction
	276
	277	Float_t curvErr=0;
	278	if(bField>0.00001){
	279	Float_t curvRadius=fMeanPtSelTrk/(0.3bField)100; //cm
	280	Float_t dRad=TMath::Sqrt(TMath::Power((gc1[0]-gc2[0]),2)+TMath::Power((gc1[1]-gc2[1]),2));
	281	Float_t aux=dRad/2.+rad1;
	282	curvErr=TMath::Sqrt(curvRadiuscurvRadius-dRaddRad/4.)-TMath::Sqrt(curvRadiuscurvRadius-auxaux); //cm
	283	}
	284
	285	Float_t sigmasq[3];
	286	sigmasq[0]=(cov[0]+curvErrcurvErr/2.)factor*factor;
	287	sigmasq[1]=(cov[3]+curvErrcurvErr/2.)factor*factor;
	288	sigmasq[2]=cov[5]factorfactor;
	289
	290	// Multiple scattering
	291	Float_t beta=1.;
	292	Float_t beta2=beta*beta;
	293	Float_t p2=fMeanPSelTrk*fMeanPSelTrk;
	294	Float_t rBP=GetPipeRadius();
	295	Float_t dBP=0.08/35.3; // 800 um of Be
	296	Float_t dL1=0.01; //approx. 1% of radiation length
	297	Float_t theta2BP=14.114.1/(beta2p21e6)TMath::Abs(dBP);
	298	Float_t theta2L1=14.114.1/(beta2p21e6)TMath::Abs(dL1);
	299	Float_t thetaBP=TMath::Sqrt(theta2BP);
	300	Float_t thetaL1=TMath::Sqrt(theta2L1);
	301	// Float_t geomfac[3];
	302	// geomfac[0]=sin(phi1)*sin(phi1);
	303	// geomfac[1]=cos(phi1)*cos(phi1);
	304	// Float_t tgth=(gc2[2]-gc1[2])/(rad2-rad1);
	305	// geomfac[2]=1+tgth*tgth;
	306	for(Int_t ico=0; ico<3;ico++){
	307	// printf("Error on coord. %d due to cov matrix+curvErr=%f\n",ico,sigmasq[ico]);
	308	// // sigmasq[ico]+=rad1rad1geomfac[ico]*theta2L1/2; // multiple scattering in layer 1
	309	// // sigmasq[ico]+=rBPrBPgeomfac[ico]*theta2BP/2; // multiple scattering in beam pipe
	310	sigmasq[ico]+=TMath::Power(rad1*TMath::Tan(thetaL1),2)/3.;
	311	sigmasq[ico]+=TMath::Power(rBP*TMath::Tan(thetaBP),2)/3.;
	312
	313	// printf("Multipl. scatt. contr %d = %f (LAY1), %f (BP)\n",ico,rad1rad1geomfac[ico]theta2L1/2,rBPrBPgeomfac[ico]theta2BP/2);
	314	// printf("Total error on coord %d = %f\n",ico,sigmasq[ico]);
	315	}
	316	Float_t wmat[9]={1.,0.,0.,0.,1.,0.,0.,0.,1.};
	317	if(sigmasq[0]!=0.) wmat[0]=1./sigmasq[0];
	318	if(sigmasq[1]!=0.) wmat[4]=1./sigmasq[1];
	319	if(sigmasq[2]!=0.) wmat[8]=1./sigmasq[2];
f5f6da22	320	new(fLines[nolines++])AliStrLine(gc1,sigmasq,wmat,gc2,kTRUE);
6a0d56b8	321
27167524	322	}
27167524	323	detTypeRec.ResetRecPoints();
70c95f95	324	}
70c95f95	325	}
70c95f95	326	}
f5f6da22	327	prpl1.Clear();
70c95f95	328	}
	329	if(nolines == 0)return -2;
	330	return nolines;
	331	}
	332
70c95f95	333	//______________________________________________________________________
05d1294c	334	Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){
70c95f95	335	// Finds the 3D vertex information using tracklets
	336	Int_t retcode = -1;
	337
ecce5370	338	Float_t xbeam=GetNominalPos()[0];
ecce5370	339	Float_t ybeam=GetNominalPos()[1];
05d1294c	340	Float_t zvert=0.;
	341	Float_t deltaR=fCoarseMaxRCut;
	342	Float_t dZmax=fZCutDiamond;
	343	if(optCuts){
27167524	344	xbeam=fVert3D.GetXv();
	345	ybeam=fVert3D.GetYv();
	346	zvert=fVert3D.GetZv();
05d1294c	347	deltaR=fMaxRCut;
	348	dZmax=fMaxZCut;
	349	}
	350
70c95f95	351	Int_t nbr=50;
05d1294c	352	Float_t rl=-fCoarseMaxRCut;
05d1294c	353	Float_t rh=fCoarseMaxRCut;
70c95f95	354	Int_t nbz=100;
	355	Float_t zl=-fZCutDiamond;
	356	Float_t zh=fZCutDiamond;
	357	Float_t binsizer=(rh-rl)/nbr;
	358	Float_t binsizez=(zh-zl)/nbz;
	359	TH3F *h3d = new TH3F("h3d","xyz distribution",nbr,rl,rh,nbr,rl,rh,nbz,zl,zh);
7340864d	360	Int_t nbrcs=25;
	361	Int_t nbzcs=50;
	362	TH3F *h3dcs = new TH3F("h3dcs","xyz distribution",nbrcs,rl,rh,nbrcs,rl,rh,nbzcs,zl,zh);
70c95f95	363
70c95f95	364	// cleanup of the TCLonesArray of tracklets (i.e. fakes are removed)
f5f6da22	365	Int_t *validate = new Int_t [fLines.GetEntriesFast()];
	366	for(Int_t i=0; i<fLines.GetEntriesFast();i++)validate[i]=0;
	367	for(Int_t i=0; i<fLines.GetEntriesFast()-1;i++){
f5f6da22	368	AliStrLine l1 = (AliStrLine)fLines.At(i);
	369	for(Int_t j=i+1;j<fLines.GetEntriesFast();j++){
	370	AliStrLine l2 = (AliStrLine)fLines.At(j);
05d1294c	371	Float_t dca=l1->GetDCA(l2);
05d1294c	372	if(dca > fDCAcut \|\| dca<0.00001) continue;
70c95f95	373	Double_t point[3];
	374	Int_t retc = l1->Cross(l2,point);
	375	if(retc<0)continue;
6a0d56b8	376	Double_t deltaZ=point[2]-zvert;
6a0d56b8	377	if(TMath::Abs(deltaZ)>dZmax)continue;
70c95f95	378	Double_t rad=TMath::Sqrt(point[0]point[0]+point[1]point[1]);
05d1294c	379	if(rad>fCoarseMaxRCut)continue;
	380	Double_t deltaX=point[0]-xbeam;
	381	Double_t deltaY=point[1]-ybeam;
05d1294c	382	Double_t raddist=TMath::Sqrt(deltaXdeltaX+deltaYdeltaY);
05d1294c	383	if(raddist>deltaR)continue;
27167524	384	validate[i]=1;
27167524	385	validate[j]=1;
70c95f95	386	h3d->Fill(point[0],point[1],point[2]);
7340864d	387	h3dcs->Fill(point[0],point[1],point[2]);
70c95f95	388	}
	389	}
	390
	391
	392
70c95f95	393	Int_t numbtracklets=0;
f5f6da22	394	for(Int_t i=0; i<fLines.GetEntriesFast();i++)if(validate[i]>=1)numbtracklets++;
e6fe443d	395	if(numbtracklets<2){delete [] validate; delete h3d; delete h3dcs; return retcode; }
70c95f95	396
f5f6da22	397	for(Int_t i=0; i<fLines.GetEntriesFast();i++){
f5f6da22	398	if(validate[i]<1)fLines.RemoveAt(i);
70c95f95	399	}
f5f6da22	400	fLines.Compress();
f5f6da22	401	AliDebug(1,Form("Number of tracklets (after compress)%d ",fLines.GetEntriesFast()));
70c95f95	402	delete [] validate;
70c95f95	403
7340864d	404	// Find peaks in histos
05d1294c	405
70c95f95	406	Double_t peak[3]={0.,0.,0.};
7340864d	407	Int_t ntrkl,ntimes;
7340864d	408	FindPeaks(h3d,peak,ntrkl,ntimes);
70c95f95	409	delete h3d;
70c95f95	410
7340864d	411	if(optCuts==0 && ntrkl<=2){
	412	ntrkl=0;
	413	ntimes=0;
	414	FindPeaks(h3dcs,peak,ntrkl,ntimes);
	415	binsizer=(rh-rl)/nbrcs;
	416	binsizez=(zh-zl)/nbzcs;
	417	if(ntrkl==1 \|\| ntimes>1){delete h3dcs; return retcode;}
	418	}
	419	delete h3dcs;
	420
70c95f95	421	// Second selection loop
7340864d	422
70c95f95	423	Float_t bs=(binsizer+binsizez)/2.;
f5f6da22	424	for(Int_t i=0; i<fLines.GetEntriesFast();i++){
	425	AliStrLine l1 = (AliStrLine)fLines.At(i);
	426	if(l1->GetDistFromPoint(peak)>2.5*bs)fLines.RemoveAt(i);
70c95f95	427	}
f5f6da22	428	fLines.Compress();
f5f6da22	429	AliDebug(1,Form("Number of tracklets (after 2nd compression) %d",fLines.GetEntriesFast()));
70c95f95	430
f5f6da22	431	if(fLines.GetEntriesFast()>1){
27167524	432	// find a first candidate for the primary vertex
f5f6da22	433	fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
70c95f95	434	// make a further selection on tracklets based on this first candidate
27167524	435	fVert3D.GetXYZ(peak);
27167524	436	AliDebug(1,Form("FIRST V candidate: %f ; %f ; %f",peak[0],peak[1],peak[2]));
f5f6da22	437	for(Int_t i=0; i<fLines.GetEntriesFast();i++){
	438	AliStrLine l1 = (AliStrLine)fLines.At(i);
	439	if(l1->GetDistFromPoint(peak)> fDCAcut)fLines.RemoveAt(i);
70c95f95	440	}
f5f6da22	441	fLines.Compress();
	442	AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines.GetEntriesFast()));
	443	if(fLines.GetEntriesFast()>1) retcode=0; // this new tracklet selection is used
27167524	444	else retcode =1; // the previous tracklet selection will be used
70c95f95	445	}
	446	else {
	447	retcode = 0;
	448	}
	449	return retcode;
	450	}
	451
6a0d56b8	452	//________________________________________________________
	453	void AliITSVertexer3D::SetMeanPPtSelTracks(Float_t fieldTesla){
	454	// Sets mean values of P and Pt based on the field
	455	if(TMath::Abs(fieldTesla-0.5)<0.01){
	456	SetMeanPSelTracks(0.885);
	457	SetMeanPtSelTracks(0.630);
	458	}else if(TMath::Abs(fieldTesla-0.4)<0.01){
	459	SetMeanPSelTracks(0.805);
	460	SetMeanPtSelTracks(0.580);
	461	}else if(TMath::Abs(fieldTesla-0.2)<0.01){
	462	SetMeanPSelTracks(0.740);
	463	SetMeanPtSelTracks(0.530);
	464	}else if(fieldTesla<0.00001){
	465	SetMeanPSelTracks(0.730);
	466	SetMeanPtSelTracks(0.510);
	467	}else{
	468	SetMeanPSelTracks();
	469	SetMeanPtSelTracks();
70c95f95	470	}
70c95f95	471	}
70c95f95	472
7340864d	473	//________________________________________________________
	474	void AliITSVertexer3D::FindPeaks(TH3F* histo, Double_t *peak, Int_t &nOfTracklets, Int_t &nOfTimes){
	475	// Finds bin with max contents in 3D histo of tracket intersections
	476	TAxis *xax = histo->GetXaxis();
	477	TAxis *yax = histo->GetYaxis();
	478	TAxis *zax = histo->GetZaxis();
	479	peak[0]=0.;
	480	peak[1]=0.;
	481	peak[2]=0.;
	482	nOfTracklets = 0;
	483	nOfTimes=0;
	484	for(Int_t i=xax->GetFirst();i<=xax->GetLast();i++){
	485	Float_t xval = xax->GetBinCenter(i);
	486	for(Int_t j=yax->GetFirst();j<=yax->GetLast();j++){
	487	Float_t yval = yax->GetBinCenter(j);
	488	for(Int_t k=zax->GetFirst();k<=zax->GetLast();k++){
	489	Float_t zval = zax->GetBinCenter(k);
	490	Int_t bc =(Int_t)histo->GetBinContent(i,j,k);
	491	if(bc>nOfTracklets){
	492	nOfTracklets = bc;
	493	peak[2] = zval;
	494	peak[1] = yval;
	495	peak[0] = xval;
	496	nOfTimes = 1;
	497	}
	498	if(bc==nOfTracklets){
	499	nOfTimes++;
	500	}
	501	}
	502	}
	503	}
	504
	505	}
70c95f95	506
	507	//________________________________________________________
	508	void AliITSVertexer3D::PrintStatus() const {
	509	// Print current status
6a0d56b8	510	printf("=======================================================\n");
	511	printf("Loose cut on Delta Phi %f\n",fCoarseDiffPhiCut);
	512	printf("Cut on tracklet DCA to Z axis %f\n",fCoarseMaxRCut);
	513	printf("Cut on tracklet DCA to beam axis %f\n",fMaxRCut);
	514	printf("Cut on diamond (Z) %f\n",fZCutDiamond);
	515	printf("Cut on DCA - tracklet to tracklet and to vertex %f\n",fDCAcut);
	516	printf(" Max Phi difference: %f\n",fDiffPhiMax);
	517	printf("=======================================================\n");
70c95f95	518	}