[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
      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.),
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
	285	AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,0,0.,0.);
	286	mv.SetTitle("Mean Vertex");
	287	mv.SetName("MeanVertex");
	288	// we have to add chi2 here
	289	AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
37add1d1	290
199775b9	291	mv.Write(mv.GetName(),TObject::kOverwrite);
	292	vtx.Write(vtx.GetName(),TObject::kOverwrite);
	293	itisOK = kTRUE;
	294	}
a3a3a28f	295	}
199775b9	296
199775b9	297	if(!itisOK){
a3a3a28f	298	AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
a3a3a28f	299	}
37add1d1	300
	301	fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
	302	fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
8e2abaa0	303	fMultH->Write(fMultH->GetName(),TObject::kOverwrite);
	304	fErrXH->Write(fErrXH->GetName(),TObject::kOverwrite);
	305	fMultHa->Write(fMultHa->GetName(),TObject::kOverwrite);
	306	fErrXHa->Write(fErrXHa->GetName(),TObject::kOverwrite);
	307	fDistH->Write(fDistH->GetName(),TObject::kOverwrite);
a1b1ad51	308	fContrH->Write(fContrH->GetName(),TObject::kOverwrite);
a1b1ad51	309	fContrHa->Write(fContrHa->GetName(),TObject::kOverwrite);
37add1d1	310	fmv.Close();
a3a3a28f	311	}
	312
	313	//______________________________________________________________________
199775b9	314	Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,UInt_t mult){
a3a3a28f	315	// Apply selection criteria to events
a3a3a28f	316	Bool_t status = kFALSE;
199775b9	317	if(!vert)return status;
8d7f6205	318	// Remove vertices reconstructed with vertexerZ
8d7f6205	319	if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
a3a3a28f	320	Int_t ncontr = vert->GetNContributors();
a3a3a28f	321	AliDebug(1,Form("Number of contributors = %d",ncontr));
199775b9	322	Double_t ex = vert->GetXRes();
8e2abaa0	323	fMultH->Fill(mult);
8e2abaa0	324	fErrXH->Fill(ex*1000.);
a1b1ad51	325	fContrH->Fill((Float_t)ncontr);
199775b9	326	if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && mult<fHighSPD0) &&
	327	((ex*1000.)<fErrXCut)) {
	328	status = kTRUE;
8e2abaa0	329	fMultHa->Fill(mult);
8e2abaa0	330	fErrXHa->Fill(ex*1000.);
a1b1ad51	331	fContrHa->Fill((Float_t)ncontr);
199775b9	332	}
a3a3a28f	333	return status;
	334	}
	335
	336	//______________________________________________________________________
	337	void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
	338	// update mean vertex
	339	Double_t currentPos[3],currentSigma[3];
	340	vert->GetXYZ(currentPos);
	341	vert->GetSigmaXYZ(currentSigma);
	342	Bool_t goon = kTRUE;
	343	for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
	344	if(!goon)return;
	345	for(Int_t i=0;i<3;i++){
762133f4	346	fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
	347	fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
	348	fAverPosSum[i]+=currentPos[i];
a3a3a28f	349	}
	350	for(Int_t i=0;i<3;i++){
	351	for(Int_t j=i;j<3;j++){
762133f4	352	fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
a3a3a28f	353	}
a3a3a28f	354	}
37add1d1	355
199775b9	356	fTotContributors+=vert->GetNContributors();
37add1d1	357	fVertexXY->Fill(currentPos[0],currentPos[1]);
	358	fVertexZ->Fill(currentPos[2]);
	359
a3a3a28f	360	fNoEventsContr++;
	361	}
	362
	363	//______________________________________________________________________
199775b9	364	Bool_t AliITSMeanVertexer::ComputeMean(Bool_t killOutliers){
	365	// compute mean vertex
	366	// called once with killOutliers = kFALSE and then again with
	367	// killOutliers = kTRUE
8e2abaa0	368
a1b1ad51	369	static Bool_t preliminary = kTRUE;
	370	static Int_t oldnumbevt = 0;
	371	if(!(preliminary \|\| killOutliers))return kTRUE; //ComputeMean is
8e2abaa0	372	// executed only once with argument kFALSE
199775b9	373	Double_t wpos[3];
199775b9	374	for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
8e2abaa0	375
a1b1ad51	376	Int_t istart = oldnumbevt;
	377	if(killOutliers)istart = 0;
	378	if(killOutliers && preliminary){
	379	preliminary = kFALSE;
	380	Reset(kTRUE,kFALSE);
8e2abaa0	381	}
a1b1ad51	382	oldnumbevt = fVertArray.GetEntries();
a1b1ad51	383	for(Int_t i=istart;i<fVertArray.GetEntries();i++){
199775b9	384	AliESDVertex* vert = NULL;
762133f4	385
199775b9	386	// second pass
	387	if(killOutliers && fAccEvents.TestBitNumber(i)){
	388	vert=(AliESDVertex*)fVertArray[i];
	389	Double_t dist=(vert->GetXv()-wpos[0])*(vert->GetXv()-wpos[0]);
	390	dist+=(vert->GetYv()-wpos[1])*(vert->GetYv()-wpos[1]);
	391	dist=sqrt(dist)*10.; // distance in mm
8e2abaa0	392	fDistH->Fill(dist);
199775b9	393	if(dist>fRCut)fAccEvents.SetBitNumber(i,kFALSE);
	394	}
	395
	396	if(!fAccEvents.TestBitNumber(i))continue;
	397	vert=(AliESDVertex*)fVertArray[i];
	398	AddToMean(vert);
	399	}
a1b1ad51	400	Bool_t bad = ((!killOutliers) && fNoEventsContr < 5) \|\| (killOutliers && fNoEventsContr <2);
	401	if(bad) {
	402	if(killOutliers){
	403	// when the control reachs this point, the preliminary evaluation of the
	404	// diamond region has to be redone. So a general reset must be issued
	405	// if this occurs, it is likely that the cuts are badly defined
	406	ResetArray();
	407	Reset(kFALSE,kTRUE);
	408	preliminary = kTRUE;
	409	oldnumbevt = 0;
	410	}
	411	return kFALSE;
	412	}
199775b9	413	Double_t nevents = fNoEventsContr;
	414	for(Int_t i=0;i<3;i++){
	415	fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i];
	416	fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
	417	fAverPos[i] = fAverPosSum[i]/nevents;
	418	}
	419	for(Int_t i=0;i<3;i++){
	420	for(Int_t j=i;j<3;j++){
	421	fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
	422	fAverPosSq[i][j] -= nevents/(nevents -1.)fAverPos[i]fAverPos[j];
	423	}
	424	}
8e2abaa0	425	if(killOutliers)ResetArray();
199775b9	426	fAverContributors = fTotContributors/nevents;
	427	return kTRUE;
	428	}