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