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

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