[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);
    fRawReader->SelectEvents(7);
  }
  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);
	135	fRawReader->SelectEvents(7);
	136	}
	137	fRunLoader = AliRunLoader::Open(fLoaderFileName.Data(),AliConfig::GetDefaultEventFolderName(),"recreate");
	138	fRunLoader->MakeTree("E");
	139	Int_t iEvent = 0;
	140	while (fRawReader->NextEvent()) {
	141	fRunLoader->SetEventNumber(iEvent);
	142	fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
	143	iEvent, iEvent);
	144	fRunLoader->MakeTree("H");
	145	fRunLoader->TreeE()->Fill();
	146	iEvent++;
	147	}
	148	fRawReader->RewindEvents();
	149	fRunLoader->SetNumberOfEventsPerFile(iEvent);
	150	fRunLoader->WriteHeader("OVERWRITE");
	151	Int_t retval = AliConfig::Instance()->AddDetector(fRunLoader->GetEventFolder(),"ITS","ITS");
	152	if(retval != 0)AliFatal("Not able to add ITS detector");
	153	AliITSLoader *loader = new AliITSLoader("ITS",fRunLoader->GetEventFolder()->GetName());
	154	fRunLoader->AddLoader(loader);
	155	fRunLoader->CdGAFile();
	156	fRunLoader->Write(0, TObject::kOverwrite);
	157	// Initialize geometry
	158
	159	AliGeomManager::LoadGeometry(fGeometryFileName.Data());
	160
	161	AliITSInitGeometry initgeom;
	162	AliITSgeom *geom = initgeom.CreateAliITSgeom();
	163	printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
	164	loader->SetITSgeom(geom);
	165	// Initialize filter values to their defaults
	166	SetFilterOnContributors();
	167	SetFilterOnTracklets();
	168	}
	169
	170	//______________________________________________________________________
	171	AliITSMeanVertexer::AliITSMeanVertexer(const AliITSMeanVertexer &vtxr) : TObject(vtxr),
	172	fLoaderFileName(vtxr.fLoaderFileName),
	173	fGeometryFileName(vtxr.fGeometryFileName),
	174	fMVFileName(vtxr.fMVFileName),
	175	fRawReader(vtxr.fRawReader),
	176	fRunLoader(vtxr.fRunLoader),
	177	fNoEventsContr(vtxr.fNoEventsContr),
	178	fTotTracklets(vtxr.fTotTracklets),
	179	fAverTracklets(vtxr.fAverTracklets),
	180	fSigmaOnAverTracks(vtxr.fSigmaOnAverTracks),
	181	fFilterOnContributors(vtxr.fFilterOnContributors),
	182	fFilterOnTracklets(vtxr.fFilterOnTracklets),
	183	fWriteVertices(vtxr.fWriteVertices)
	184	{
	185	// Copy constructor
	186	// Copies are not allowed. The method is protected to avoid misuse.
	187	AliFatal("Copy constructor not allowed\n");
	188
189	}
190
191	//______________________________________________________________________
192	AliITSMeanVertexer& AliITSMeanVertexer::operator=(const AliITSMeanVertexer& /* vtxr */ ){
193	// Assignment operator
194	AliError("Assignment operator not allowed\n");
195	return *this;
196	}
197
198	//______________________________________________________________________
199	AliITSMeanVertexer::~AliITSMeanVertexer() {
200	// Destructor
201	delete fRawReader;
202	delete fRunLoader;
203
204	}
205
206	//______________________________________________________________________
207	void AliITSMeanVertexer::Reconstruct(){
208	// Performs SPD local reconstruction
209	AliITSLoader* loader = static_cast<AliITSLoader*>(fRunLoader->GetLoader("ITSLoader"));
210	if (!loader) {
211	AliFatal("ITS loader not found");
212	return;
213	}
214	AliITSDetTypeRec* rec = new AliITSDetTypeRec();
215	rec->SetITSgeom(loader->GetITSgeom());
216	rec->SetDefaults();
217
218	rec->SetDefaultClusterFindersV2(kTRUE);
219
220	Int_t nEvents = fRunLoader->GetNumberOfEvents();
221	fRawReader->RewindEvents();
222	for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
223	fRawReader->NextEvent();
224	fRunLoader->GetEvent(iEvent);
225	AliDebug(1,Form(">>>>>>> Processing event number: %d",iEvent));
226	loader->LoadRecPoints("update");
227	loader->CleanRecPoints();
228	loader->MakeRecPointsContainer();
229	TTree *tR = loader->TreeR();
230	if(!tR){
231	AliFatal("Tree R pointer not found - Abort \n");
232	break;
233	}
234	rec->DigitsToRecPoints(fRawReader,tR,"SPD");
235	rec->ResetRecPoints();
236	rec->ResetClusters();
237	loader->WriteRecPoints("OVERWRITE");
238	loader->UnloadRecPoints();
239	}
240
241	}
242
243	//______________________________________________________________________
244	void AliITSMeanVertexer::DoVertices(){
245	// Loop on all events and compute 3D vertices
246	AliITSLoader* loader = static_cast<AliITSLoader*>(fRunLoader->GetLoader("ITSLoader"));
308c2f7c	247	AliITSVertexer3D *dovert = new AliITSVertexer3D();
a3a3a28f	248	AliESDVertex *vert = 0;
	249	Int_t nevents = fRunLoader->TreeE()->GetEntries();
	250	for(Int_t i=0; i<nevents; i++){
	251	fRunLoader->GetEvent(i);
308c2f7c	252	TTree* cltree = loader->TreeR();
308c2f7c	253	vert = dovert->FindVertexForCurrentEvent(cltree);
a3a3a28f	254	AliMultiplicity *mult = dovert->GetMultiplicity();
	255	if(Filter(vert,mult)){
	256	AddToMean(vert);
	257	if(fWriteVertices){
	258	loader->PostVertex(vert);
	259	loader->WriteVertices();
	260	}
	261	}
	262	else {
	263	if(vert)delete vert;
	264	}
	265	}
	266	if(ComputeMean()){
	267	Double_t cov[6];
	268	cov[0] = fAverPosSq[0][0]; // variance x
	269	cov[1] = fAverPosSq[0][1]; // cov xy
	270	cov[2] = fAverPosSq[1][1]; // variance y
	271	cov[3] = fAverPosSq[0][2]; // cov xz
	272	cov[4] = fAverPosSq[1][2]; // cov yz
	273	cov[5] = fAverPosSq[2][2]; // variance z
	274	AliMeanVertex mv(fWeighPos,fWeighSig,cov,nevents,fTotTracklets,fAverTracklets,fSigmaOnAverTracks);
	275	mv.SetTitle("Mean Vertex");
	276	mv.SetName("Meanvertex");
	277	AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets()));
	278	AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks()));
	279	TFile fmv(fMVFileName.Data(),"recreate");
	280	mv.Write();
	281	fmv.Close();
	282	}
	283	else {
	284	AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
	285	}
	286	delete dovert;
	287	}
	288
	289	//______________________________________________________________________
	290	Bool_t AliITSMeanVertexer::Filter(AliESDVertex vert,AliMultiplicity mult){
	291	// Apply selection criteria to events
	292	Bool_t status = kFALSE;
	293	if(!vert \|\| !mult)return status;
	294	Int_t ncontr = vert->GetNContributors();
	295	Int_t ntracklets = mult->GetNumberOfTracklets();
	296	AliDebug(1,Form("Number of contributors = %d",ncontr));
	297	AliDebug(1,Form("Number of tracklets = %d",ntracklets));
	298	if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE;
	299	fAverTracklets += ntracklets;
	300	fSigmaOnAverTracks += ntracklets*ntracklets;
	301	return status;
	302	}
	303
	304	//______________________________________________________________________
	305	void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
	306	// update mean vertex
	307	Double_t currentPos[3],currentSigma[3];
	308	vert->GetXYZ(currentPos);
	309	vert->GetSigmaXYZ(currentSigma);
	310	Bool_t goon = kTRUE;
	311	for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
	312	if(!goon)return;
	313	for(Int_t i=0;i<3;i++){
	314	fWeighPos[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
	315	fWeighSig[i]+=1./currentSigma[i]/currentSigma[i];
	316	fAverPos[i]+=currentPos[i];
	317	}
318	for(Int_t i=0;i<3;i++){
319	for(Int_t j=i;j<3;j++){
320	fAverPosSq[i][j] += currentPos[i] * currentPos[j];
321	}
322	}
323	fNoEventsContr++;
324	}
325
326	//______________________________________________________________________
327	Bool_t AliITSMeanVertexer::ComputeMean(){
328	// compute mean vertex
329	if(fNoEventsContr < 2) return kFALSE;
330	Double_t nevents = fNoEventsContr;
331	for(Int_t i=0;i<3;i++){
332	fWeighPos[i] /= fWeighSig[i];
333	fWeighSig[i] = 1./TMath::Sqrt(fWeighSig[i]);
334	fAverPos[i] /= nevents;
335	}
336	for(Int_t i=0;i<3;i++){
337	for(Int_t j=i;j<3;j++){
338	fAverPosSq[i][j] /= (nevents -1.);
339	fAverPosSq[i][j] -= nevents/(nevents -1.)fAverPos[i]fAverPos[j];
340	}
341	}
342	fTotTracklets = fAverTracklets; // total number of tracklets used
343	fAverTracklets /= nevents;
344	fSigmaOnAverTracks /= (nevents - 1);
345	Double_t tmp = nevents/(nevents -1.)fAverTrackletsfAverTracklets;
346	fSigmaOnAverTracks -= tmp;
347	fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks);
348	return kTRUE;
349	}