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

#include <TFile.h>
#include <TH1.h>
#include <TH2.h>
#include "AliGeomManager.h"
#include "AliITSDetTypeRec.h"
#include "AliITSInitGeometry.h"
#include "AliITSMeanVertexer.h"
#include "AliITSRecPointContainer.h"
#include "AliITSLoader.h"
#include "AliLog.h"
#include "AliRawReader.h"
#include "AliRawReaderDate.h"
#include "AliRawReaderRoot.h"
#include "AliRunLoader.h"
#include "AliITSVertexer3D.h"
#include "AliITSVertexer3DTapan.h"
#include "AliESDVertex.h"
#include "AliMeanVertex.h"

const Int_t AliITSMeanVertexer::fgkMaxNumOfEvents = 10000;
ClassImp(AliITSMeanVertexer)

///////////////////////////////////////////////////////////////////////
//                                                                   //
// Class to compute vertex position using SPD local reconstruction   //
// An average vertex position using all the events                   //
// is built and saved                                                //
// Individual vertices can be optionally stored                      //
// Origin: M.Masera  (masera@to.infn.it)                             //
// Usage:
// AliITSMeanVertexer mv("RawDataFileName");
// mv.SetGeometryFileName("FileWithGeometry.root"); // def. geometry.root
// ....  optional setters ....
// mv.Reconstruct();  // SPD local reconstruction
// mv.DoVertices(); 
// Resulting AliMeanVertex object is stored in a file named fMVFileName
///////////////////////////////////////////////////////////////////////

/* $Id$ */

//______________________________________________________________________
AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(),
fDetTypeRec(NULL),
fVertexXY(NULL),
fVertexZ(NULL),
fNoEventsContr(0),
fTotContributors(0.),
fAverContributors(0.), 
fFilterOnContributors(0),
fMode(mode),
fVertexer(NULL),
fAccEvents(fgkMaxNumOfEvents), 
fVertArray("AliESDVertex",fgkMaxNumOfEvents),
fClu0(NULL),
fIndex(0),
fErrXCut(0.),
fRCut(0.),
fLowSPD0(0),
fHighSPD0(0)
{
  // Default Constructor
  Reset();

  // Histograms initialization
  const Float_t xLimit = 5.0, yLimit = 5.0, zLimit = 50.0;
  const Float_t xDelta = 0.02, yDelta = 0.02, zDelta = 0.2;
  fVertexXY = new TH2F("VertexXY","Vertex Diamond (Y vs X)",
		       2*(Int_t)(xLimit/xDelta),-xLimit,xLimit,
		       2*(Int_t)(yLimit/yDelta),-yLimit,yLimit);
  fVertexXY->SetXTitle("X , cm");
  fVertexXY->SetYTitle("Y , cm");
  fVertexXY->SetOption("colz");
  fVertexZ  = new TH1F("VertexZ"," Longitudinal Vertex Profile",
		       2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
  fVertexZ->SetXTitle("Z , cm");
  fClu0 = new UInt_t [fgkMaxNumOfEvents];
  for(Int_t i=0;i<fgkMaxNumOfEvents;i++)fClu0[i]=0;
  fAccEvents.ResetAllBits(kTRUE);
}

//______________________________________________________________________
void AliITSMeanVertexer::Reset(){
  // Reset averages 
  for(Int_t i=0;i<3;i++){
    fWeighPosSum[i] = 0.;
    fWeighSigSum[i] = 0.;
    fAverPosSum[i] = 0.;
    fWeighPos[i] = 0.;
    fWeighSig[i] = 0.;
    fAverPos[i] = 0.;
    for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
    for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
    fNoEventsContr=0;
    fTotContributors = 0.;
    if(fVertexXY)fVertexXY->Reset();
    if(fVertexZ)fVertexZ->Reset();
  }
}
//______________________________________________________________________
Bool_t AliITSMeanVertexer::Init() {
  // Initialize filters
  // Initialize geometry
  // Initialize ITS classes
 
  AliGeomManager::LoadGeometry();
  if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;

  AliITSInitGeometry initgeom;
  AliITSgeom *geom = initgeom.CreateAliITSgeom();
  if (!geom) return kFALSE;
  printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());

  fDetTypeRec = new AliITSDetTypeRec();
  fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
  fDetTypeRec->SetITSgeom(geom);
  fDetTypeRec->SetDefaults();
  fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);

  // Initialize filter values to their defaults
  SetFilterOnContributors();
  SetCutOnErrX();
  SetCutOnR();
  SetCutOnCls();

  // Instatiate vertexer
  if (!fMode) {
    fVertexer = new AliITSVertexer3DTapan(1000);
  }
  else {
    fVertexer = new AliITSVertexer3D();
    fVertexer->SetDetTypeRec(fDetTypeRec);
    AliITSVertexer3D* alias = (AliITSVertexer3D*)fVertexer;
    alias->SetWideFiducialRegion(40.,0.5);
    alias->SetNarrowFiducialRegion(0.5,0.5);
    alias->SetDeltaPhiCuts(0.5,0.025);
    alias->SetDCACut(0.1);
    alias->SetPileupAlgo(3);
    fVertexer->SetComputeMultiplicity(kFALSE);
  }
  return kTRUE;
}

//______________________________________________________________________
AliITSMeanVertexer::~AliITSMeanVertexer() {
  // Destructor
  delete fDetTypeRec;
  delete fVertexXY;
  delete fVertexZ;
  delete fVertexer;
}

//______________________________________________________________________
Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader){
  // Performs SPD local reconstruction
  // and vertex finding
  // returns true in case a vertex is found

  // Run SPD cluster finder
  AliITSRecPointContainer::Instance()->PrepareToRead();
  TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
  fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");

  Bool_t vtxOK = kFALSE;
  UInt_t nl1=0;
  AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
  if (!fMode) {
    if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
  }
  else {
    AliITSRecPointContainer* rpcont= AliITSRecPointContainer::Instance();
    nl1=rpcont->GetNClustersInLayerFast(1);
    if(Filter(vtx,nl1)) vtxOK = kTRUE;
  }
  delete clustersTree;
  if (vtxOK && fMode){
    new(fVertArray[fIndex])AliESDVertex(*vtx);
    fClu0[fIndex]=nl1;
    fIndex++;
  }
  if (vtx) delete vtx;

  return vtxOK;
}

//______________________________________________________________________
void AliITSMeanVertexer::WriteVertices(const char *filename){
  // Compute mean vertex and
  // store it along with the histograms
  // in a file
  
  TFile fmv(filename,"update");
  Bool_t itisOK = kFALSE;
  if(ComputeMean(kFALSE)){
    if(ComputeMean(kTRUE)){
      Double_t cov[6];
      cov[0] =  fAverPosSq[0][0];  // variance x
      cov[1] =  fAverPosSq[0][1];  // cov xy
      cov[2] =  fAverPosSq[1][1];  // variance y
      cov[3] =  fAverPosSq[0][2];  // cov xz
      cov[4] =  fAverPosSq[1][2];  // cov yz
      cov[5] =  fAverPosSq[2][2];  // variance z
      AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
      mv.SetTitle("Mean Vertex");
      mv.SetName("MeanVertex");
      // we have to add chi2 here
      AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");

      mv.Write(mv.GetName(),TObject::kOverwrite);
      vtx.Write(vtx.GetName(),TObject::kOverwrite);
      itisOK = kTRUE;
    }
  }

  if(!itisOK){
    AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
  }

  fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
  fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
  fmv.Close();
}

//______________________________________________________________________
Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,UInt_t mult){
  // Apply selection criteria to events
  Bool_t status = kFALSE;
  if(!vert)return status;
  // Remove vertices reconstructed with vertexerZ
  if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
  Int_t ncontr = vert->GetNContributors();
  AliDebug(1,Form("Number of contributors = %d",ncontr));
  Double_t ex = vert->GetXRes();
  if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && mult<fHighSPD0) && 
     ((ex*1000.)<fErrXCut)) {
    status = kTRUE;
    fAccEvents.SetBitNumber(fIndex);
  }
  return status;
}

//______________________________________________________________________
void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
  // update mean vertex
  Double_t currentPos[3],currentSigma[3];
  vert->GetXYZ(currentPos);
  vert->GetSigmaXYZ(currentSigma);
  Bool_t goon = kTRUE;
  for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
  if(!goon)return;
  for(Int_t i=0;i<3;i++){
    fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
    fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
    fAverPosSum[i]+=currentPos[i];
  }
  for(Int_t i=0;i<3;i++){
    for(Int_t j=i;j<3;j++){
      fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
    }
  }

  fTotContributors+=vert->GetNContributors();
  fVertexXY->Fill(currentPos[0],currentPos[1]);
  fVertexZ->Fill(currentPos[2]);

  fNoEventsContr++;
}

//______________________________________________________________________
 Bool_t AliITSMeanVertexer::ComputeMean(Bool_t killOutliers){
   // compute mean vertex
   // called once with killOutliers = kFALSE and then again with
   // killOutliers = kTRUE
   Double_t wpos[3];
   for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
   Reset();
   for(Int_t i=0;i<fVertArray.GetEntries();i++){
     AliESDVertex* vert = NULL;

     // second pass
     if(killOutliers && fAccEvents.TestBitNumber(i)){
       vert=(AliESDVertex*)fVertArray[i];
       Double_t dist=(vert->GetXv()-wpos[0])*(vert->GetXv()-wpos[0]);
       dist+=(vert->GetYv()-wpos[1])*(vert->GetYv()-wpos[1]);
       dist=sqrt(dist)*10.;    // distance in mm
       if(dist>fRCut)fAccEvents.SetBitNumber(i,kFALSE);
     }
     
     if(!fAccEvents.TestBitNumber(i))continue;
     vert=(AliESDVertex*)fVertArray[i];
     AddToMean(vert);
   }
   if(fNoEventsContr < 2) return kFALSE;
   Double_t nevents = fNoEventsContr;
   for(Int_t i=0;i<3;i++){
     fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i]; 
     fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
     fAverPos[i] = fAverPosSum[i]/nevents;
   } 
   for(Int_t i=0;i<3;i++){
     for(Int_t j=i;j<3;j++){
       fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
       fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
     }
   } 
   
   fAverContributors = fTotContributors/nevents;
   return kTRUE;
 }
Commit	Line	Data
	1	#include <TFile.h>
	2	#include <TH1.h>
	3	#include <TH2.h>
	4	#include "AliGeomManager.h"
	5	#include "AliITSDetTypeRec.h"
	6	#include "AliITSInitGeometry.h"
	7	#include "AliITSMeanVertexer.h"
	8	#include "AliITSRecPointContainer.h"
	9	#include "AliITSLoader.h"
	10	#include "AliLog.h"
	11	#include "AliRawReader.h"
	12	#include "AliRawReaderDate.h"
	13	#include "AliRawReaderRoot.h"
	14	#include "AliRunLoader.h"
	15	#include "AliITSVertexer3D.h"
	16	#include "AliITSVertexer3DTapan.h"
	17	#include "AliESDVertex.h"
	18	#include "AliMeanVertex.h"
	19
	20	const Int_t AliITSMeanVertexer::fgkMaxNumOfEvents = 10000;
	21	ClassImp(AliITSMeanVertexer)
	22
	23	///////////////////////////////////////////////////////////////////////
	24	// //
	25	// Class to compute vertex position using SPD local reconstruction //
	26	// An average vertex position using all the events //
	27	// is built and saved //
	28	// Individual vertices can be optionally stored //
	29	// Origin: M.Masera (masera@to.infn.it) //
	30	// Usage:
	31	// AliITSMeanVertexer mv("RawDataFileName");
	32	// mv.SetGeometryFileName("FileWithGeometry.root"); // def. geometry.root
	33	// .... optional setters ....
	34	// mv.Reconstruct(); // SPD local reconstruction
	35	// mv.DoVertices();
	36	// Resulting AliMeanVertex object is stored in a file named fMVFileName
	37	///////////////////////////////////////////////////////////////////////
	38
	39	/* $Id$ */
	40
	41	//______________________________________________________________________
	42	AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(),
	43	fDetTypeRec(NULL),
	44	fVertexXY(NULL),
	45	fVertexZ(NULL),
	46	fNoEventsContr(0),
	47	fTotContributors(0.),
	48	fAverContributors(0.),
	49	fFilterOnContributors(0),
	50	fMode(mode),
	51	fVertexer(NULL),
	52	fAccEvents(fgkMaxNumOfEvents),
	53	fVertArray("AliESDVertex",fgkMaxNumOfEvents),
	54	fClu0(NULL),
	55	fIndex(0),
	56	fErrXCut(0.),
	57	fRCut(0.),
	58	fLowSPD0(0),
	59	fHighSPD0(0)
	60	{
	61	// Default Constructor
	62	Reset();
	63
	64	// Histograms initialization
	65	const Float_t xLimit = 5.0, yLimit = 5.0, zLimit = 50.0;
	66	const Float_t xDelta = 0.02, yDelta = 0.02, zDelta = 0.2;
	67	fVertexXY = new TH2F("VertexXY","Vertex Diamond (Y vs X)",
	68	2*(Int_t)(xLimit/xDelta),-xLimit,xLimit,
	69	2*(Int_t)(yLimit/yDelta),-yLimit,yLimit);
	70	fVertexXY->SetXTitle("X , cm");
	71	fVertexXY->SetYTitle("Y , cm");
	72	fVertexXY->SetOption("colz");
	73	fVertexZ = new TH1F("VertexZ"," Longitudinal Vertex Profile",
	74	2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
	75	fVertexZ->SetXTitle("Z , cm");
	76	fClu0 = new UInt_t [fgkMaxNumOfEvents];
	77	for(Int_t i=0;i<fgkMaxNumOfEvents;i++)fClu0[i]=0;
	78	fAccEvents.ResetAllBits(kTRUE);
	79	}
	80
	81	//______________________________________________________________________
	82	void AliITSMeanVertexer::Reset(){
	83	// Reset averages
	84	for(Int_t i=0;i<3;i++){
	85	fWeighPosSum[i] = 0.;
	86	fWeighSigSum[i] = 0.;
	87	fAverPosSum[i] = 0.;
	88	fWeighPos[i] = 0.;
	89	fWeighSig[i] = 0.;
	90	fAverPos[i] = 0.;
	91	for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
	92	for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
	93	fNoEventsContr=0;
	94	fTotContributors = 0.;
	95	if(fVertexXY)fVertexXY->Reset();
	96	if(fVertexZ)fVertexZ->Reset();
	97	}
	98	}
	99	//______________________________________________________________________
	100	Bool_t AliITSMeanVertexer::Init() {
	101	// Initialize filters
	102	// Initialize geometry
	103	// Initialize ITS classes
	104
	105	AliGeomManager::LoadGeometry();
	106	if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
	107
	108	AliITSInitGeometry initgeom;
	109	AliITSgeom *geom = initgeom.CreateAliITSgeom();
	110	if (!geom) return kFALSE;
	111	printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
	112
	113	fDetTypeRec = new AliITSDetTypeRec();
	114	fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
	115	fDetTypeRec->SetITSgeom(geom);
	116	fDetTypeRec->SetDefaults();
	117	fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
	118
	119	// Initialize filter values to their defaults
	120	SetFilterOnContributors();
	121	SetCutOnErrX();
	122	SetCutOnR();
	123	SetCutOnCls();
	124
	125	// Instatiate vertexer
	126	if (!fMode) {
	127	fVertexer = new AliITSVertexer3DTapan(1000);
	128	}
	129	else {
	130	fVertexer = new AliITSVertexer3D();
	131	fVertexer->SetDetTypeRec(fDetTypeRec);
	132	AliITSVertexer3D* alias = (AliITSVertexer3D*)fVertexer;
	133	alias->SetWideFiducialRegion(40.,0.5);
	134	alias->SetNarrowFiducialRegion(0.5,0.5);
	135	alias->SetDeltaPhiCuts(0.5,0.025);
	136	alias->SetDCACut(0.1);
	137	alias->SetPileupAlgo(3);
	138	fVertexer->SetComputeMultiplicity(kFALSE);
	139	}
	140	return kTRUE;
	141	}
	142
	143	//______________________________________________________________________
	144	AliITSMeanVertexer::~AliITSMeanVertexer() {
	145	// Destructor
	146	delete fDetTypeRec;
	147	delete fVertexXY;
	148	delete fVertexZ;
	149	delete fVertexer;
	150	}
	151
	152	//______________________________________________________________________
	153	Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader){
	154	// Performs SPD local reconstruction
	155	// and vertex finding
	156	// returns true in case a vertex is found
	157
	158	// Run SPD cluster finder
	159	AliITSRecPointContainer::Instance()->PrepareToRead();
	160	TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
	161	fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
	162
	163	Bool_t vtxOK = kFALSE;
	164	UInt_t nl1=0;
	165	AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
	166	if (!fMode) {
	167	if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
	168	}
	169	else {
	170	AliITSRecPointContainer* rpcont= AliITSRecPointContainer::Instance();
	171	nl1=rpcont->GetNClustersInLayerFast(1);
	172	if(Filter(vtx,nl1)) vtxOK = kTRUE;
	173	}
	174	delete clustersTree;
	175	if (vtxOK && fMode){
	176	new(fVertArray[fIndex])AliESDVertex(*vtx);
	177	fClu0[fIndex]=nl1;
	178	fIndex++;
	179	}
	180	if (vtx) delete vtx;
	181
	182	return vtxOK;
	183	}
	184
	185	//______________________________________________________________________
	186	void AliITSMeanVertexer::WriteVertices(const char *filename){
	187	// Compute mean vertex and
	188	// store it along with the histograms
	189	// in a file
	190
	191	TFile fmv(filename,"update");
	192	Bool_t itisOK = kFALSE;
	193	if(ComputeMean(kFALSE)){
	194	if(ComputeMean(kTRUE)){
	195	Double_t cov[6];
	196	cov[0] = fAverPosSq[0][0]; // variance x
	197	cov[1] = fAverPosSq[0][1]; // cov xy
	198	cov[2] = fAverPosSq[1][1]; // variance y
	199	cov[3] = fAverPosSq[0][2]; // cov xz
	200	cov[4] = fAverPosSq[1][2]; // cov yz
	201	cov[5] = fAverPosSq[2][2]; // variance z
	202	AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
	203	mv.SetTitle("Mean Vertex");
	204	mv.SetName("MeanVertex");
	205	// we have to add chi2 here
	206	AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
	207
	208	mv.Write(mv.GetName(),TObject::kOverwrite);
	209	vtx.Write(vtx.GetName(),TObject::kOverwrite);
	210	itisOK = kTRUE;
	211	}
	212	}
	213
	214	if(!itisOK){
	215	AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
	216	}
	217
	218	fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
	219	fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
	220	fmv.Close();
	221	}
	222
	223	//______________________________________________________________________
	224	Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,UInt_t mult){
	225	// Apply selection criteria to events
	226	Bool_t status = kFALSE;
	227	if(!vert)return status;
	228	// Remove vertices reconstructed with vertexerZ
	229	if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
	230	Int_t ncontr = vert->GetNContributors();
	231	AliDebug(1,Form("Number of contributors = %d",ncontr));
	232	Double_t ex = vert->GetXRes();
	233	if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && mult<fHighSPD0) &&
	234	((ex*1000.)<fErrXCut)) {
	235	status = kTRUE;
	236	fAccEvents.SetBitNumber(fIndex);
	237	}
	238	return status;
	239	}
	240
	241	//______________________________________________________________________
	242	void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
	243	// update mean vertex
	244	Double_t currentPos[3],currentSigma[3];
	245	vert->GetXYZ(currentPos);
	246	vert->GetSigmaXYZ(currentSigma);
	247	Bool_t goon = kTRUE;
	248	for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
	249	if(!goon)return;
	250	for(Int_t i=0;i<3;i++){
	251	fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
	252	fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
	253	fAverPosSum[i]+=currentPos[i];
	254	}
	255	for(Int_t i=0;i<3;i++){
	256	for(Int_t j=i;j<3;j++){
	257	fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
	258	}
	259	}
	260
	261	fTotContributors+=vert->GetNContributors();
	262	fVertexXY->Fill(currentPos[0],currentPos[1]);
	263	fVertexZ->Fill(currentPos[2]);
	264
	265	fNoEventsContr++;
	266	}
	267
	268	//______________________________________________________________________
	269	Bool_t AliITSMeanVertexer::ComputeMean(Bool_t killOutliers){
	270	// compute mean vertex
	271	// called once with killOutliers = kFALSE and then again with
	272	// killOutliers = kTRUE
	273	Double_t wpos[3];
	274	for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
	275	Reset();
	276	for(Int_t i=0;i<fVertArray.GetEntries();i++){
	277	AliESDVertex* vert = NULL;
	278
	279	// second pass
	280	if(killOutliers && fAccEvents.TestBitNumber(i)){
	281	vert=(AliESDVertex*)fVertArray[i];
	282	Double_t dist=(vert->GetXv()-wpos[0])*(vert->GetXv()-wpos[0]);
	283	dist+=(vert->GetYv()-wpos[1])*(vert->GetYv()-wpos[1]);
	284	dist=sqrt(dist)*10.; // distance in mm
	285	if(dist>fRCut)fAccEvents.SetBitNumber(i,kFALSE);
	286	}
	287
	288	if(!fAccEvents.TestBitNumber(i))continue;
	289	vert=(AliESDVertex*)fVertArray[i];
	290	AddToMean(vert);
	291	}
	292	if(fNoEventsContr < 2) return kFALSE;
	293	Double_t nevents = fNoEventsContr;
	294	for(Int_t i=0;i<3;i++){
	295	fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i];
	296	fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
	297	fAverPos[i] = fAverPosSum[i]/nevents;
	298	}
	299	for(Int_t i=0;i<3;i++){
	300	for(Int_t j=i;j<3;j++){
	301	fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
	302	fAverPosSq[i][j] -= nevents/(nevents -1.)fAverPos[i]fAverPos[j];
	303	}
	304	}
	305
	306	fAverContributors = fTotContributors/nevents;
	307	return kTRUE;
	308	}