[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"
//#include "AliCodeTimer.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:                                                            //
// Class used by the ITSSPDVertexDiamondda.cxx detector algorithm    //
///////////////////////////////////////////////////////////////////////

/* $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),
fMultH(NULL),
fErrXH(NULL),
fMultHa(NULL),
fErrXHa(NULL),
fDistH(NULL)
{
  // 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");
  fErrXH = new TH1F("errorX","Error X - before cuts",100,0.,99.);
  fMultH = new TH1F("multh","mult on layer 1 before cuts",2400,0.,7200.);
  fErrXHa = new TH1F("errorXa","Error X - after Filter",100,0.,99.);
  fErrXHa->SetLineColor(kRed);
  fMultHa = new TH1F("multha","mult on layer 1 - after Filter",2400,0.,7200.);
  fMultHa->SetLineColor(kRed);
  fDistH = new TH1F("disth","distance (mm)",100,0.,30.);
  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 fMultH;
  delete fErrXH;
  delete fMultHa;
  delete fErrXHa;
  delete fDistH;
  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
  static Int_t evcount = -1;
  if(evcount <0){
    evcount++;
    AliInfo(Form("Low and high cuts on SPD L0 clusters %d , %d \n",fLowSPD0,fHighSPD0));
    AliInfo(Form("Reconstruct: cut on errX %f \n",fErrXCut));
  }
//  AliCodeTimerAuto("",0);
  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;
    /*    if(vtx){
      if(vtxOK){
	printf("The vertex is OK\n");
      }
      else {
	printf("The vertex is NOT OK\n");
      }
      vtx->PrintStatus();
    }
    else {
      printf("The vertex was not reconstructed\n");
      }  */
  }
  delete clustersTree;
  if (vtxOK && fMode){
    new(fVertArray[fIndex])AliESDVertex(*vtx);
    fClu0[fIndex]=nl1;
    fAccEvents.SetBitNumber(fIndex);
    fIndex++;
    if(fIndex>=fgkMaxNumOfEvents){
      if(ComputeMean(kFALSE)){
	if(ComputeMean(kTRUE))AliInfo("Mean vertex computed");
      }
    }
  }
  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);
  fMultH->Write(fMultH->GetName(),TObject::kOverwrite);
  fErrXH->Write(fErrXH->GetName(),TObject::kOverwrite);
  fMultHa->Write(fMultHa->GetName(),TObject::kOverwrite);
  fErrXHa->Write(fErrXHa->GetName(),TObject::kOverwrite);
  fDistH->Write(fDistH->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();
  fMultH->Fill(mult);
  fErrXH->Fill(ex*1000.);
  if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && mult<fHighSPD0) && 
     ((ex*1000.)<fErrXCut)) {
    status = kTRUE;
    fMultHa->Fill(mult);
    fErrXHa->Fill(ex*1000.);
  }
  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

   static Bool_t firstime = kTRUE;
   if(fNoEventsContr>2 && !killOutliers)return kTRUE;  //ComputeMean is
                             // executed only once with argument kFALSE
   Double_t wpos[3];
   for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];

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