[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
a3a3a28f	1	#include <TFile.h>
37add1d1	2	#include <TH1.h>
37add1d1	3	#include <TH2.h>
a3a3a28f	4	#include "AliGeomManager.h"
a3a3a28f	5	#include "AliITSDetTypeRec.h"
	6	#include "AliITSInitGeometry.h"
	7	#include "AliITSMeanVertexer.h"
d5eae813	8	#include "AliITSRecPointContainer.h"
a3a3a28f	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"
37add1d1	16	#include "AliITSVertexer3DTapan.h"
a3a3a28f	17	#include "AliESDVertex.h"
a3a3a28f	18	#include "AliMeanVertex.h"
a3a3a28f	19
199775b9	20	const Int_t AliITSMeanVertexer::fgkMaxNumOfEvents = 10000;
a3a3a28f	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	//______________________________________________________________________
d5eae813	42	AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(),
37add1d1	43	fDetTypeRec(NULL),
	44	fVertexXY(NULL),
	45	fVertexZ(NULL),
a3a3a28f	46	fNoEventsContr(0),
199775b9	47	fTotContributors(0.),
199775b9	48	fAverContributors(0.),
a3a3a28f	49	fFilterOnContributors(0),
d5eae813	50	fMode(mode),
199775b9	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)
a3a3a28f	60	{
a3a3a28f	61	// Default Constructor
199775b9	62	Reset();
a3a3a28f	63
37add1d1	64	// Histograms initialization
8d7f6205	65	const Float_t xLimit = 5.0, yLimit = 5.0, zLimit = 50.0;
37add1d1	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);
63845149	70	fVertexXY->SetXTitle("X , cm");
	71	fVertexXY->SetYTitle("Y , cm");
	72	fVertexXY->SetOption("colz");
37add1d1	73	fVertexZ = new TH1F("VertexZ"," Longitudinal Vertex Profile",
37add1d1	74	2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
63845149	75	fVertexZ->SetXTitle("Z , cm");
199775b9	76	fClu0 = new UInt_t [fgkMaxNumOfEvents];
	77	for(Int_t i=0;i<fgkMaxNumOfEvents;i++)fClu0[i]=0;
	78	fAccEvents.ResetAllBits(kTRUE);
a3a3a28f	79	}
a3a3a28f	80
199775b9	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	}
a3a3a28f	99	//______________________________________________________________________
37add1d1	100	Bool_t AliITSMeanVertexer::Init() {
37add1d1	101	// Initialize filters
a3a3a28f	102	// Initialize geometry
37add1d1	103	// Initialize ITS classes
a3a3a28f	104
37add1d1	105	AliGeomManager::LoadGeometry();
37add1d1	106	if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
a3a3a28f	107
	108	AliITSInitGeometry initgeom;
	109	AliITSgeom *geom = initgeom.CreateAliITSgeom();
37add1d1	110	if (!geom) return kFALSE;
a3a3a28f	111	printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
37add1d1	112
cf0b1cea	113	fDetTypeRec = new AliITSDetTypeRec();
762133f4	114	fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
37add1d1	115	fDetTypeRec->SetITSgeom(geom);
	116	fDetTypeRec->SetDefaults();
	117	fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
	118
a3a3a28f	119	// Initialize filter values to their defaults
a3a3a28f	120	SetFilterOnContributors();
199775b9	121	SetCutOnErrX();
	122	SetCutOnR();
	123	SetCutOnCls();
a3a3a28f	124
d5eae813	125	// Instatiate vertexer
	126	if (!fMode) {
	127	fVertexer = new AliITSVertexer3DTapan(1000);
	128	}
	129	else {
	130	fVertexer = new AliITSVertexer3D();
	131	fVertexer->SetDetTypeRec(fDetTypeRec);
199775b9	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);
d5eae813	139	}
37add1d1	140	return kTRUE;
a3a3a28f	141	}
	142
	143	//______________________________________________________________________
	144	AliITSMeanVertexer::~AliITSMeanVertexer() {
	145	// Destructor
37add1d1	146	delete fDetTypeRec;
	147	delete fVertexXY;
	148	delete fVertexZ;
d5eae813	149	delete fVertexer;
a3a3a28f	150	}
	151
	152	//______________________________________________________________________
d5eae813	153	Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader){
a3a3a28f	154	// Performs SPD local reconstruction
37add1d1	155	// and vertex finding
	156	// returns true in case a vertex is found
	157
	158	// Run SPD cluster finder
d5eae813	159	AliITSRecPointContainer::Instance()->PrepareToRead();
37add1d1	160	TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
	161	fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
	162
	163	Bool_t vtxOK = kFALSE;
199775b9	164	UInt_t nl1=0;
d5eae813	165	AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
d5eae813	166	if (!fMode) {
37add1d1	167	if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
a3a3a28f	168	}
37add1d1	169	else {
199775b9	170	AliITSRecPointContainer* rpcont= AliITSRecPointContainer::Instance();
	171	nl1=rpcont->GetNClustersInLayerFast(1);
	172	if(Filter(vtx,nl1)) vtxOK = kTRUE;
a3a3a28f	173	}
37add1d1	174	delete clustersTree;
199775b9	175	if (vtxOK && fMode){
	176	new(fVertArray[fIndex])AliESDVertex(*vtx);
	177	fClu0[fIndex]=nl1;
	178	fIndex++;
	179	}
37add1d1	180	if (vtx) delete vtx;
a3a3a28f	181
37add1d1	182	return vtxOK;
a3a3a28f	183	}
	184
	185	//______________________________________________________________________
37add1d1	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");
199775b9	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");
37add1d1	207
199775b9	208	mv.Write(mv.GetName(),TObject::kOverwrite);
	209	vtx.Write(vtx.GetName(),TObject::kOverwrite);
	210	itisOK = kTRUE;
	211	}
a3a3a28f	212	}
199775b9	213
199775b9	214	if(!itisOK){
a3a3a28f	215	AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
a3a3a28f	216	}
37add1d1	217
	218	fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
	219	fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
	220	fmv.Close();
a3a3a28f	221	}
	222
	223	//______________________________________________________________________
199775b9	224	Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,UInt_t mult){
a3a3a28f	225	// Apply selection criteria to events
a3a3a28f	226	Bool_t status = kFALSE;
199775b9	227	if(!vert)return status;
8d7f6205	228	// Remove vertices reconstructed with vertexerZ
8d7f6205	229	if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
a3a3a28f	230	Int_t ncontr = vert->GetNContributors();
a3a3a28f	231	AliDebug(1,Form("Number of contributors = %d",ncontr));
199775b9	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	}
a3a3a28f	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++){
762133f4	251	fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
	252	fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
	253	fAverPosSum[i]+=currentPos[i];
a3a3a28f	254	}
	255	for(Int_t i=0;i<3;i++){
	256	for(Int_t j=i;j<3;j++){
762133f4	257	fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
a3a3a28f	258	}
a3a3a28f	259	}
37add1d1	260
199775b9	261	fTotContributors+=vert->GetNContributors();
37add1d1	262	fVertexXY->Fill(currentPos[0],currentPos[1]);
	263	fVertexZ->Fill(currentPos[2]);
	264
a3a3a28f	265	fNoEventsContr++;
	266	}
	267
	268	//______________________________________________________________________
199775b9	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;
762133f4	278
199775b9	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	}