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