[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.),
fZCutDiamond(0.),
fLowSPD0(0),
fHighSPD0(0),
fMultH(NULL),
fErrXH(NULL),
fMultHa(NULL),
fErrXHa(NULL),
fDistH(NULL),
fContrH(NULL),
fContrHa(NULL)
{
  // Default Constructor
  SetZFiducialRegion();   //  sets fZCutDiamond to its default value
  Reset(kFALSE,kFALSE);

  // 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.);
  fContrH = new TH1F("contrib","Number of contributors - before cuts",300,0.5,300.5);
  fContrHa = new TH1F("contriba","Number of contributors - after cuts",300,0.5,300.5);
  fContrHa->SetLineColor(kRed);

  fClu0 = new UInt_t [fgkMaxNumOfEvents];
  for(Int_t i=0;i<fgkMaxNumOfEvents;i++)fClu0[i]=0;
  fAccEvents.ResetAllBits(kTRUE);
}

//______________________________________________________________________
void AliITSMeanVertexer::Reset(Bool_t redefine2D,Bool_t complete){
  // Reset averages 
  // Both arguments are kFALSE when the method is called by the constructor
  // redefine2D is TRUE when the method is called by ComputeMean at the second
  //            pass.
  // redefine2D is FALSE and complete is TRUE when the method is called
  //            after a failure cof ComputeMean(kTRUE)

  static Int_t counter = 0;
  if(fVertexXY && redefine2D){
    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);
  }
  else if(fVertexXY && !redefine2D){
    fVertexXY->Reset();
  }
  if(fVertexZ){
    fVertexZ->Reset();
    fDistH->Reset();
    if(complete){
      fErrXH->Reset();
      fMultH->Reset();
      fErrXHa->Reset(); 
      fMultHa->Reset();
      fContrH->Reset();
      fContrHa->Reset();
    }
  }
  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.;
  }
}
//______________________________________________________________________
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(fZCutDiamond,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 fContrH;
  delete fContrHa;
  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
      // We use standard average and not weighed averag now
      // AliMeanVertex is apparently not taken by the preprocessor; only
      // the AliESDVertex object is retrieved
      //      AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
      AliMeanVertex mv(fAverPos,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");
      AliESDVertex vtx(fAverPos,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);
  fContrH->Write(fContrH->GetName(),TObject::kOverwrite);
  fContrHa->Write(fContrHa->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.);
  fContrH->Fill((Float_t)ncontr);
  if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && mult<fHighSPD0) && 
     ((ex*1000.)<fErrXCut)) {
    status = kTRUE;
    fMultHa->Fill(mult);
    fErrXHa->Fill(ex*1000.);
    fContrHa->Fill((Float_t)ncontr);
  }
  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 preliminary = kTRUE;
   static Int_t oldnumbevt = 0;
   if(!(preliminary || 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];

   Int_t istart = oldnumbevt;
   if(killOutliers)istart = 0;
   if(killOutliers && preliminary){
     preliminary = kFALSE;
     Reset(kTRUE,kFALSE);
   }
   oldnumbevt = fVertArray.GetEntries();
   for(Int_t i=istart;i<fVertArray.GetEntries();i++){
     AliESDVertex* vert = NULL;

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