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

#include <TFile.h>
#include "AliGeomManager.h"
#include "AliHeader.h"
#include "AliITSDetTypeRec.h"
#include "AliITSInitGeometry.h"
#include "AliITSMeanVertexer.h"
#include "AliITSLoader.h"
#include "AliLog.h"
#include "AliRawReader.h"
#include "AliRawReaderDate.h"
#include "AliRawReaderRoot.h"
#include "AliRunLoader.h"
#include "AliITSVertexer3D.h"
#include "AliESDVertex.h"
#include "AliMeanVertex.h"
#include "AliMultiplicity.h"

const TString AliITSMeanVertexer::fgkMVFileNameDefault = "MeanVertex.root";


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():TObject(),
fLoaderFileName(),
fGeometryFileName(),
fMVFileName(),
fRawReader(),
fRunLoader(),
fNoEventsContr(0),
fTotTracklets(0.),
fAverTracklets(0.),
fSigmaOnAverTracks(0.), 
fFilterOnContributors(0),
fFilterOnTracklets(0),
fWriteVertices(kFALSE)
{
  // Default Constructor
  for(Int_t i=0;i<3;i++){
    fWeighPos[i] = 0.;
    fWeighSig[i] = 0.;
    fAverPos[i] = 0.;
    for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
  }
  
}

//______________________________________________________________________
AliITSMeanVertexer::AliITSMeanVertexer(TString &filename):TObject(),
fLoaderFileName(),
fGeometryFileName(),
fMVFileName(fgkMVFileNameDefault),
fRawReader(),
fRunLoader(),
fNoEventsContr(0),
fTotTracklets(0.),
fAverTracklets(0.),
fSigmaOnAverTracks(0.), 
fFilterOnContributors(0),
fFilterOnTracklets(0),
fWriteVertices(kTRUE)
{
  // Standard constructor

  for(Int_t i=0;i<3;i++){
    fWeighPos[i] = 0.;
    fWeighSig[i] = 0.;
    fAverPos[i] = 0.;
    for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
  }
  SetLoaderFileName();
  SetGeometryFileName();

  Init(filename);
}
//______________________________________________________________________
AliITSMeanVertexer::AliITSMeanVertexer(TString &filename, 
                                       TString &loaderfilename, 
				       TString &geometryfilename):TObject(),
fLoaderFileName(),
fGeometryFileName(),
fMVFileName(fgkMVFileNameDefault),
fRawReader(),
fRunLoader(),
fNoEventsContr(0), 
fTotTracklets(0.),
fAverTracklets(0.),
fSigmaOnAverTracks(0.), 
fFilterOnContributors(0),
fFilterOnTracklets(0),
fWriteVertices(kTRUE)
{
  // Standard constructor with explicit geometry file name assignment
  for(Int_t i=0;i<3;i++){
    fWeighPos[i] = 0.;
    fWeighSig[i] = 0.;
    fAverPos[i] = 0.;
    for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
  }
  SetLoaderFileName(loaderfilename);
  SetGeometryFileName(geometryfilename);
  Init(filename);
}

//______________________________________________________________________
void AliITSMeanVertexer::Init(TString &filename){
  // Initialization part common to different constructors
  if(filename.IsNull()){
    AliFatal("Please, provide a valid file name for raw data file\n");
  }
  // if file name ends with root a raw reader ROOT is assumed
  if(filename.EndsWith(".root")){
    fRawReader = new AliRawReaderRoot(filename);
  }
  else {  // DATE raw reader is assumed
    fRawReader = new AliRawReaderDate(filename);
  }
  fRunLoader = AliRunLoader::Open(fLoaderFileName.Data(),AliConfig::GetDefaultEventFolderName(),"recreate");
  fRunLoader->MakeTree("E");
  Int_t iEvent = 0;
  while (fRawReader->NextEvent()) {
    fRunLoader->SetEventNumber(iEvent);
    fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(), 
				   iEvent, iEvent);
    fRunLoader->MakeTree("H");
    fRunLoader->TreeE()->Fill();
    iEvent++;
  }
  fRawReader->RewindEvents();
  fRunLoader->SetNumberOfEventsPerFile(iEvent);
  fRunLoader->WriteHeader("OVERWRITE");
 Int_t retval = AliConfig::Instance()->AddDetector(fRunLoader->GetEventFolder(),"ITS","ITS");
 if(retval != 0)AliFatal("Not able to add ITS detector");
  AliITSLoader *loader = new AliITSLoader("ITS",fRunLoader->GetEventFolder()->GetName());
  fRunLoader->AddLoader(loader);
  fRunLoader->CdGAFile();
  fRunLoader->Write(0, TObject::kOverwrite);
  // Initialize geometry
 
  AliGeomManager::LoadGeometry(fGeometryFileName.Data());

  AliITSInitGeometry initgeom;
  AliITSgeom *geom = initgeom.CreateAliITSgeom();
  printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
  loader->SetITSgeom(geom);
  // Initialize filter values to their defaults
  SetFilterOnContributors();
  SetFilterOnTracklets();
}

//______________________________________________________________________
AliITSMeanVertexer::AliITSMeanVertexer(const AliITSMeanVertexer &vtxr) : TObject(vtxr),
fLoaderFileName(vtxr.fLoaderFileName),
fGeometryFileName(vtxr.fGeometryFileName),
fMVFileName(vtxr.fMVFileName),
fRawReader(vtxr.fRawReader),
fRunLoader(vtxr.fRunLoader),
fNoEventsContr(vtxr.fNoEventsContr),
fTotTracklets(vtxr.fTotTracklets),
fAverTracklets(vtxr.fAverTracklets),
fSigmaOnAverTracks(vtxr.fSigmaOnAverTracks),  
fFilterOnContributors(vtxr.fFilterOnContributors),
fFilterOnTracklets(vtxr.fFilterOnTracklets),
fWriteVertices(vtxr.fWriteVertices)
{
  // Copy constructor
  // Copies are not allowed. The method is protected to avoid misuse.
  AliFatal("Copy constructor not allowed\n");

}

//______________________________________________________________________
AliITSMeanVertexer& AliITSMeanVertexer::operator=(const AliITSMeanVertexer&  /* vtxr */ ){
  // Assignment operator
  AliError("Assignment operator not allowed\n");
  return *this;
}

//______________________________________________________________________
AliITSMeanVertexer::~AliITSMeanVertexer() {
  // Destructor
  delete fRawReader;
  delete fRunLoader;

}

//______________________________________________________________________
void AliITSMeanVertexer::Reconstruct(){
  // Performs SPD local reconstruction
  AliITSLoader* loader = static_cast<AliITSLoader*>(fRunLoader->GetLoader("ITSLoader"));
  if (!loader) {
    AliFatal("ITS loader not found");
    return;
  }
  AliITSDetTypeRec* rec = new AliITSDetTypeRec();
  rec->SetITSgeom(loader->GetITSgeom());
  rec->SetDefaults();

  rec->SetDefaultClusterFindersV2(kTRUE);

  Int_t nEvents = fRunLoader->GetNumberOfEvents();
  fRawReader->RewindEvents();
  for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
    fRawReader->NextEvent();
    fRunLoader->GetEvent(iEvent);
    AliDebug(1,Form(">>>>>>>   Processing event number: %d",iEvent));
    loader->LoadRecPoints("update");
    loader->CleanRecPoints();
    loader->MakeRecPointsContainer();
    TTree *tR = loader->TreeR();
    if(!tR){
      AliFatal("Tree R pointer not found - Abort \n");
      break;
    }
    rec->DigitsToRecPoints(fRawReader,tR,"SPD");
    rec->ResetRecPoints();
    rec->ResetClusters();    
    loader->WriteRecPoints("OVERWRITE");
    loader->UnloadRecPoints();
  }

}

//______________________________________________________________________
void AliITSMeanVertexer::DoVertices(){
  // Loop on all events and compute 3D vertices
  AliITSLoader* loader = static_cast<AliITSLoader*>(fRunLoader->GetLoader("ITSLoader"));
  AliITSVertexer3D *dovert = new AliITSVertexer3D();
  AliESDVertex *vert = 0;
  Int_t nevents = fRunLoader->TreeE()->GetEntries();
  for(Int_t i=0; i<nevents; i++){
    fRunLoader->GetEvent(i);
    TTree* cltree = loader->TreeR();
    vert = dovert->FindVertexForCurrentEvent(cltree);
    AliMultiplicity *mult = dovert->GetMultiplicity();
    if(Filter(vert,mult)){
      AddToMean(vert); 
      if(fWriteVertices){
	loader->PostVertex(vert);
	loader->WriteVertices();
      }
    }
    else {
      if(vert)delete vert;
    }
  }
  if(ComputeMean()){
    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,nevents,fTotTracklets,fAverTracklets,fSigmaOnAverTracks);
    mv.SetTitle("Mean Vertex");
    mv.SetName("Meanvertex");
    AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets()));
    AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks()));
    TFile fmv(fMVFileName.Data(),"recreate");
    mv.Write();
    fmv.Close();
  }
  else {
    AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
  }
  delete dovert;
}

//______________________________________________________________________
Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,AliMultiplicity *mult){
  // Apply selection criteria to events
  Bool_t status = kFALSE;
  if(!vert || !mult)return status;
  Int_t ncontr = vert->GetNContributors();
  Int_t ntracklets = mult->GetNumberOfTracklets();
  AliDebug(1,Form("Number of contributors = %d",ncontr));
  AliDebug(1,Form("Number of tracklets = %d",ntracklets));
  if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE;
  fAverTracklets += ntracklets;
  fSigmaOnAverTracks += ntracklets*ntracklets;
  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++){
    fWeighPos[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
    fWeighSig[i]+=1./currentSigma[i]/currentSigma[i];
    fAverPos[i]+=currentPos[i];
  }
  for(Int_t i=0;i<3;i++){
    for(Int_t j=i;j<3;j++){
      fAverPosSq[i][j] += currentPos[i] * currentPos[j];
    }
  }
  fNoEventsContr++;
}

//______________________________________________________________________
Bool_t AliITSMeanVertexer::ComputeMean(){
  // compute mean vertex
  if(fNoEventsContr < 2) return kFALSE;
  Double_t nevents = fNoEventsContr;
  for(Int_t i=0;i<3;i++){
    fWeighPos[i] /= fWeighSig[i]; 
    fWeighSig[i] = 1./TMath::Sqrt(fWeighSig[i]);
    fAverPos[i] /= nevents;
  } 
  for(Int_t i=0;i<3;i++){
    for(Int_t j=i;j<3;j++){
      fAverPosSq[i][j] /= (nevents -1.);
      fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
    }
  } 
  fTotTracklets = fAverTracklets;  //  total number of tracklets used 
  fAverTracklets /= nevents;
  fSigmaOnAverTracks /= (nevents - 1);
  Double_t tmp = nevents/(nevents -1.)*fAverTracklets*fAverTracklets;
  fSigmaOnAverTracks -= tmp;
  fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks);
  return kTRUE;
}
Commit	Line	Data
a3a3a28f	1	#include <TFile.h>
	2	#include "AliGeomManager.h"
	3	#include "AliHeader.h"
	4	#include "AliITSDetTypeRec.h"
	5	#include "AliITSInitGeometry.h"
	6	#include "AliITSMeanVertexer.h"
	7	#include "AliITSLoader.h"
	8	#include "AliLog.h"
	9	#include "AliRawReader.h"
	10	#include "AliRawReaderDate.h"
	11	#include "AliRawReaderRoot.h"
	12	#include "AliRunLoader.h"
	13	#include "AliITSVertexer3D.h"
	14	#include "AliESDVertex.h"
	15	#include "AliMeanVertex.h"
	16	#include "AliMultiplicity.h"
	17
	18	const TString AliITSMeanVertexer::fgkMVFileNameDefault = "MeanVertex.root";
	19
	20
	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():TObject(),
	43	fLoaderFileName(),
	44	fGeometryFileName(),
	45	fMVFileName(),
	46	fRawReader(),
	47	fRunLoader(),
	48	fNoEventsContr(0),
	49	fTotTracklets(0.),
	50	fAverTracklets(0.),
	51	fSigmaOnAverTracks(0.),
	52	fFilterOnContributors(0),
	53	fFilterOnTracklets(0),
	54	fWriteVertices(kFALSE)
	55	{
	56	// Default Constructor
	57	for(Int_t i=0;i<3;i++){
	58	fWeighPos[i] = 0.;
	59	fWeighSig[i] = 0.;
	60	fAverPos[i] = 0.;
	61	for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
	62	}
	63
	64	}
65
66	//______________________________________________________________________
308c2f7c	67	AliITSMeanVertexer::AliITSMeanVertexer(TString &filename):TObject(),
a3a3a28f	68	fLoaderFileName(),
	69	fGeometryFileName(),
	70	fMVFileName(fgkMVFileNameDefault),
	71	fRawReader(),
	72	fRunLoader(),
	73	fNoEventsContr(0),
	74	fTotTracklets(0.),
	75	fAverTracklets(0.),
	76	fSigmaOnAverTracks(0.),
	77	fFilterOnContributors(0),
	78	fFilterOnTracklets(0),
	79	fWriteVertices(kTRUE)
	80	{
	81	// Standard constructor
	82
	83	for(Int_t i=0;i<3;i++){
	84	fWeighPos[i] = 0.;
	85	fWeighSig[i] = 0.;
	86	fAverPos[i] = 0.;
	87	for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
	88	}
	89	SetLoaderFileName();
	90	SetGeometryFileName();
	91
	92	Init(filename);
	93	}
	94	//______________________________________________________________________
308c2f7c	95	AliITSMeanVertexer::AliITSMeanVertexer(TString &filename,
	96	TString &loaderfilename,
	97	TString &geometryfilename):TObject(),
a3a3a28f	98	fLoaderFileName(),
	99	fGeometryFileName(),
	100	fMVFileName(fgkMVFileNameDefault),
	101	fRawReader(),
	102	fRunLoader(),
	103	fNoEventsContr(0),
	104	fTotTracklets(0.),
	105	fAverTracklets(0.),
	106	fSigmaOnAverTracks(0.),
	107	fFilterOnContributors(0),
	108	fFilterOnTracklets(0),
	109	fWriteVertices(kTRUE)
	110	{
	111	// Standard constructor with explicit geometry file name assignment
	112	for(Int_t i=0;i<3;i++){
	113	fWeighPos[i] = 0.;
	114	fWeighSig[i] = 0.;
	115	fAverPos[i] = 0.;
	116	for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
	117	}
	118	SetLoaderFileName(loaderfilename);
	119	SetGeometryFileName(geometryfilename);
	120	Init(filename);
	121	}
	122
	123	//______________________________________________________________________
308c2f7c	124	void AliITSMeanVertexer::Init(TString &filename){
a3a3a28f	125	// Initialization part common to different constructors
	126	if(filename.IsNull()){
	127	AliFatal("Please, provide a valid file name for raw data file\n");
	128	}
	129	// if file name ends with root a raw reader ROOT is assumed
	130	if(filename.EndsWith(".root")){
	131	fRawReader = new AliRawReaderRoot(filename);
	132	}
	133	else { // DATE raw reader is assumed
	134	fRawReader = new AliRawReaderDate(filename);
a3a3a28f	135	}
	136	fRunLoader = AliRunLoader::Open(fLoaderFileName.Data(),AliConfig::GetDefaultEventFolderName(),"recreate");
	137	fRunLoader->MakeTree("E");
	138	Int_t iEvent = 0;
	139	while (fRawReader->NextEvent()) {
	140	fRunLoader->SetEventNumber(iEvent);
	141	fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
	142	iEvent, iEvent);
	143	fRunLoader->MakeTree("H");
	144	fRunLoader->TreeE()->Fill();
	145	iEvent++;
	146	}
	147	fRawReader->RewindEvents();
	148	fRunLoader->SetNumberOfEventsPerFile(iEvent);
	149	fRunLoader->WriteHeader("OVERWRITE");
	150	Int_t retval = AliConfig::Instance()->AddDetector(fRunLoader->GetEventFolder(),"ITS","ITS");
	151	if(retval != 0)AliFatal("Not able to add ITS detector");
	152	AliITSLoader *loader = new AliITSLoader("ITS",fRunLoader->GetEventFolder()->GetName());
	153	fRunLoader->AddLoader(loader);
	154	fRunLoader->CdGAFile();
	155	fRunLoader->Write(0, TObject::kOverwrite);
	156	// Initialize geometry
	157
	158	AliGeomManager::LoadGeometry(fGeometryFileName.Data());
	159
	160	AliITSInitGeometry initgeom;
	161	AliITSgeom *geom = initgeom.CreateAliITSgeom();
	162	printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
	163	loader->SetITSgeom(geom);
	164	// Initialize filter values to their defaults
	165	SetFilterOnContributors();
	166	SetFilterOnTracklets();
	167	}
	168
	169	//______________________________________________________________________
	170	AliITSMeanVertexer::AliITSMeanVertexer(const AliITSMeanVertexer &vtxr) : TObject(vtxr),
	171	fLoaderFileName(vtxr.fLoaderFileName),
	172	fGeometryFileName(vtxr.fGeometryFileName),
	173	fMVFileName(vtxr.fMVFileName),
	174	fRawReader(vtxr.fRawReader),
	175	fRunLoader(vtxr.fRunLoader),
	176	fNoEventsContr(vtxr.fNoEventsContr),
	177	fTotTracklets(vtxr.fTotTracklets),
	178	fAverTracklets(vtxr.fAverTracklets),
	179	fSigmaOnAverTracks(vtxr.fSigmaOnAverTracks),
	180	fFilterOnContributors(vtxr.fFilterOnContributors),
	181	fFilterOnTracklets(vtxr.fFilterOnTracklets),
	182	fWriteVertices(vtxr.fWriteVertices)
	183	{
	184	// Copy constructor
	185	// Copies are not allowed. The method is protected to avoid misuse.
	186	AliFatal("Copy constructor not allowed\n");
	187
	188	}
	189
	190	//______________________________________________________________________
	191	AliITSMeanVertexer& AliITSMeanVertexer::operator=(const AliITSMeanVertexer& /* vtxr */ ){
	192	// Assignment operator
	193	AliError("Assignment operator not allowed\n");
	194	return *this;
	195	}
	196
	197	//______________________________________________________________________
	198	AliITSMeanVertexer::~AliITSMeanVertexer() {
199	// Destructor
200	delete fRawReader;
201	delete fRunLoader;
202
203	}
204
205	//______________________________________________________________________
206	void AliITSMeanVertexer::Reconstruct(){
207	// Performs SPD local reconstruction
208	AliITSLoader* loader = static_cast<AliITSLoader*>(fRunLoader->GetLoader("ITSLoader"));
209	if (!loader) {
210	AliFatal("ITS loader not found");
211	return;
212	}
213	AliITSDetTypeRec* rec = new AliITSDetTypeRec();
214	rec->SetITSgeom(loader->GetITSgeom());
215	rec->SetDefaults();
216
217	rec->SetDefaultClusterFindersV2(kTRUE);
218
219	Int_t nEvents = fRunLoader->GetNumberOfEvents();
220	fRawReader->RewindEvents();
221	for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
222	fRawReader->NextEvent();
223	fRunLoader->GetEvent(iEvent);
224	AliDebug(1,Form(">>>>>>> Processing event number: %d",iEvent));
225	loader->LoadRecPoints("update");
226	loader->CleanRecPoints();
227	loader->MakeRecPointsContainer();
228	TTree *tR = loader->TreeR();
229	if(!tR){
230	AliFatal("Tree R pointer not found - Abort \n");
231	break;
232	}
233	rec->DigitsToRecPoints(fRawReader,tR,"SPD");
234	rec->ResetRecPoints();
235	rec->ResetClusters();
236	loader->WriteRecPoints("OVERWRITE");
237	loader->UnloadRecPoints();
238	}
239
240	}
241
242	//______________________________________________________________________
243	void AliITSMeanVertexer::DoVertices(){
244	// Loop on all events and compute 3D vertices
245	AliITSLoader* loader = static_cast<AliITSLoader*>(fRunLoader->GetLoader("ITSLoader"));
308c2f7c	246	AliITSVertexer3D *dovert = new AliITSVertexer3D();
a3a3a28f	247	AliESDVertex *vert = 0;
	248	Int_t nevents = fRunLoader->TreeE()->GetEntries();
	249	for(Int_t i=0; i<nevents; i++){
	250	fRunLoader->GetEvent(i);
308c2f7c	251	TTree* cltree = loader->TreeR();
308c2f7c	252	vert = dovert->FindVertexForCurrentEvent(cltree);
a3a3a28f	253	AliMultiplicity *mult = dovert->GetMultiplicity();
	254	if(Filter(vert,mult)){
	255	AddToMean(vert);
	256	if(fWriteVertices){
	257	loader->PostVertex(vert);
	258	loader->WriteVertices();
	259	}
	260	}
	261	else {
	262	if(vert)delete vert;
	263	}
	264	}
	265	if(ComputeMean()){
	266	Double_t cov[6];
	267	cov[0] = fAverPosSq[0][0]; // variance x
	268	cov[1] = fAverPosSq[0][1]; // cov xy
	269	cov[2] = fAverPosSq[1][1]; // variance y
	270	cov[3] = fAverPosSq[0][2]; // cov xz
	271	cov[4] = fAverPosSq[1][2]; // cov yz
	272	cov[5] = fAverPosSq[2][2]; // variance z
	273	AliMeanVertex mv(fWeighPos,fWeighSig,cov,nevents,fTotTracklets,fAverTracklets,fSigmaOnAverTracks);
	274	mv.SetTitle("Mean Vertex");
	275	mv.SetName("Meanvertex");
	276	AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets()));
	277	AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks()));
	278	TFile fmv(fMVFileName.Data(),"recreate");
	279	mv.Write();
	280	fmv.Close();
	281	}
	282	else {
	283	AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
	284	}
	285	delete dovert;
	286	}
	287
	288	//______________________________________________________________________
	289	Bool_t AliITSMeanVertexer::Filter(AliESDVertex vert,AliMultiplicity mult){
	290	// Apply selection criteria to events
	291	Bool_t status = kFALSE;
	292	if(!vert \|\| !mult)return status;
	293	Int_t ncontr = vert->GetNContributors();
	294	Int_t ntracklets = mult->GetNumberOfTracklets();
	295	AliDebug(1,Form("Number of contributors = %d",ncontr));
	296	AliDebug(1,Form("Number of tracklets = %d",ntracklets));
	297	if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE;
	298	fAverTracklets += ntracklets;
	299	fSigmaOnAverTracks += ntracklets*ntracklets;
	300	return status;
	301	}
	302
	303	//______________________________________________________________________
	304	void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
	305	// update mean vertex
	306	Double_t currentPos[3],currentSigma[3];
	307	vert->GetXYZ(currentPos);
	308	vert->GetSigmaXYZ(currentSigma);
	309	Bool_t goon = kTRUE;
	310	for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
	311	if(!goon)return;
	312	for(Int_t i=0;i<3;i++){
	313	fWeighPos[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
	314	fWeighSig[i]+=1./currentSigma[i]/currentSigma[i];
	315	fAverPos[i]+=currentPos[i];
	316	}
317	for(Int_t i=0;i<3;i++){
318	for(Int_t j=i;j<3;j++){
319	fAverPosSq[i][j] += currentPos[i] * currentPos[j];
320	}
321	}
322	fNoEventsContr++;
323	}
324
325	//______________________________________________________________________
326	Bool_t AliITSMeanVertexer::ComputeMean(){
327	// compute mean vertex
328	if(fNoEventsContr < 2) return kFALSE;
329	Double_t nevents = fNoEventsContr;
330	for(Int_t i=0;i<3;i++){
331	fWeighPos[i] /= fWeighSig[i];
332	fWeighSig[i] = 1./TMath::Sqrt(fWeighSig[i]);
333	fAverPos[i] /= nevents;
334	}
335	for(Int_t i=0;i<3;i++){
336	for(Int_t j=i;j<3;j++){
337	fAverPosSq[i][j] /= (nevents -1.);
338	fAverPosSq[i][j] -= nevents/(nevents -1.)fAverPos[i]fAverPos[j];
339	}
340	}
341	fTotTracklets = fAverTracklets; // total number of tracklets used
342	fAverTracklets /= nevents;
343	fSigmaOnAverTracks /= (nevents - 1);
344	Double_t tmp = nevents/(nevents -1.)fAverTrackletsfAverTracklets;
345	fSigmaOnAverTracks -= tmp;
346	fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks);
347	return kTRUE;
348	}