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