[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->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
	1	#include <TFile.h>
	2	#include <TH1.h>
	3	#include <TH2.h>
	4	#include <TTree.h>
	5	#include "AliGeomManager.h"
	6	#include "AliITSDetTypeRec.h"
	7	#include "AliITSInitGeometry.h"
	8	#include "AliITSMeanVertexer.h"
	9	#include "AliITSRecPointContainer.h"
	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"
	17	#include "AliITSVertexer3DTapan.h"
	18	#include "AliESDVertex.h"
	19	#include "AliMeanVertex.h"
	20	//#include "AliCodeTimer.h"
	21
	22	const Int_t AliITSMeanVertexer::fgkMaxNumOfEvents = 10000;
	23	ClassImp(AliITSMeanVertexer)
	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) //
	31	// Usage: //
	32	// Class used by the ITSSPDVertexDiamondda.cxx detector algorithm //
	33	///////////////////////////////////////////////////////////////////////
	34	/* $Id$ */
	35	//______________________________________________________________________
	36	AliITSMeanVertexer::AliITSMeanVertexer(Bool_t mode):TObject(),
	37	fDetTypeRec(NULL),
	38	fVertexXY(NULL),
	39	fVertexZ(NULL),
	40	fNoEventsContr(0),
	41	fTotContributors(0.),
	42	fAverContributors(0.),
	43	fFilterOnContributors(0),
	44	fMode(mode),
	45	fVertexer(NULL),
	46	fAccEvents(fgkMaxNumOfEvents),
	47	fVertArray("AliESDVertex",fgkMaxNumOfEvents),
	48	fClu0(NULL),
	49	fIndex(0),
	50	fErrXCut(0.),
	51	fRCut(0.),
	52	fZCutDiamond(0.),
	53	fLowSPD0(0),
	54	fHighSPD0(0),
	55	fMultH(NULL),
	56	fErrXH(NULL),
	57	fMultHa(NULL),
	58	fErrXHa(NULL),
	59	fDistH(NULL),
	60	fContrH(NULL),
	61	fContrHa(NULL)
	62	{
	63	// Default Constructor
	64	SetZFiducialRegion(); // sets fZCutDiamond to its default value
	65	Reset(kFALSE,kFALSE);
	66
	67	// Histograms initialization
	68	const Float_t xLimit = 0.6, yLimit = 0.6, zLimit = 50.0;
	69	const Float_t xDelta = 0.003, yDelta = 0.003, zDelta = 0.2;
	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);
	73	fVertexXY->SetXTitle("X , cm");
	74	fVertexXY->SetYTitle("Y , cm");
	75	fVertexXY->SetOption("colz");
	76	fVertexZ = new TH1F("VertexZ"," Longitudinal Vertex Profile",
	77	2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
	78	fVertexZ->SetXTitle("Z , cm");
	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.);
	82	fErrXHa->SetLineColor(kRed);
	83	fMultHa = new TH1F("multha","mult on layer 1 - after Filter",2400,0.,7200.);
	84	fMultHa->SetLineColor(kRed);
	85	fDistH = new TH1F("disth","distance (mm)",100,0.,30.);
	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
	90	fClu0 = new UInt_t [fgkMaxNumOfEvents];
	91	for(Int_t i=0;i<fgkMaxNumOfEvents;i++)fClu0[i]=0;
	92	fAccEvents.ResetAllBits(kTRUE);
	93	}
	94
	95	//______________________________________________________________________
	96	void AliITSMeanVertexer::Reset(Bool_t redefine2D,Bool_t complete){
	97	// Reset averages
	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
	104	static Int_t counter = 0;
	105	if(fVertexXY && redefine2D){
	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	}
	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	}
	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.;
	148	}
	149	}
	150	//______________________________________________________________________
	151	Bool_t AliITSMeanVertexer::Init() {
	152	// Initialize filters
	153	// Initialize geometry
	154	// Initialize ITS classes
	155
	156	AliGeomManager::LoadGeometry();
	157	if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
	158
	159	AliITSInitGeometry initgeom;
	160	AliITSgeom *geom = initgeom.CreateAliITSgeom();
	161	if (!geom) return kFALSE;
	162	printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
	163
	164	fDetTypeRec = new AliITSDetTypeRec();
	165	fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
	166	fDetTypeRec->SetITSgeom(geom);
	167	fDetTypeRec->SetDefaults();
	168	fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
	169
	170	// Initialize filter values to their defaults
	171	SetFilterOnContributors();
	172	SetCutOnErrX();
	173	SetCutOnR();
	174	SetCutOnCls();
	175
	176	// Instatiate vertexer
	177	if (!fMode) {
	178	fVertexer = new AliITSVertexer3DTapan(1000);
	179	}
	180	else {
	181	fVertexer = new AliITSVertexer3D();
	182	fVertexer->SetDetTypeRec(fDetTypeRec);
	183	AliITSVertexer3D* alias = (AliITSVertexer3D*)fVertexer;
	184	alias->SetWideFiducialRegion(fZCutDiamond,0.5);
	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);
	190	}
	191	return kTRUE;
	192	}
	193
	194	//______________________________________________________________________
	195	AliITSMeanVertexer::~AliITSMeanVertexer() {
	196	// Destructor
	197	delete fDetTypeRec;
	198	delete fVertexXY;
	199	delete fVertexZ;
	200	delete fMultH;
	201	delete fErrXH;
	202	delete fMultHa;
	203	delete fErrXHa;
	204	delete fDistH;
	205	delete fContrH;
	206	delete fContrHa;
	207	delete fVertexer;
	208	}
	209
	210	//______________________________________________________________________
	211	Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader){
	212	// Performs SPD local reconstruction
	213	// and vertex finding
	214	// returns true in case a vertex is found
	215
	216	// Run SPD cluster finder
	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	}
	223	// AliCodeTimerAuto("",0);
	224	AliITSRecPointContainer::Instance()->PrepareToRead();
	225	TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
	226	fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
	227
	228	Bool_t vtxOK = kFALSE;
	229	UInt_t nl1=0;
	230	AliESDVertex *vtx = fVertexer->FindVertexForCurrentEvent(clustersTree);
	231	if (!fMode) {
	232	if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
	233	}
	234	else {
	235	AliITSRecPointContainer* rpcont= AliITSRecPointContainer::Instance();
	236	nl1=rpcont->GetNClustersInLayerFast(1);
	237	if(Filter(vtx,nl1)) vtxOK = kTRUE;
	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	} */
	250	}
	251	delete clustersTree;
	252	if (vtxOK && fMode){
	253	new(fVertArray[fIndex])AliESDVertex(*vtx);
	254	fClu0[fIndex]=nl1;
	255	fAccEvents.SetBitNumber(fIndex);
	256	fIndex++;
	257	if(fIndex>=fgkMaxNumOfEvents){
	258	if(ComputeMean(kFALSE)){
	259	if(ComputeMean(kTRUE))AliInfo("Mean vertex computed");
	260	}
	261	}
	262	}
	263	if (vtx) delete vtx;
	264
	265	return vtxOK;
	266	}
	267
	268	//______________________________________________________________________
	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");
	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	// 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.);
	290	mv.SetTitle("Mean Vertex");
	291	mv.SetName("MeanVertex");
	292	// we have to add chi2 here
	293	// AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
	294	AliESDVertex vtx(fAverPos,cov,0,TMath::Nint(fAverContributors),"MeanVertexPos");
	295
	296	mv.Write(mv.GetName(),TObject::kOverwrite);
	297	vtx.Write(vtx.GetName(),TObject::kOverwrite);
	298	itisOK = kTRUE;
	299	}
	300	}
	301
	302	if(!itisOK){
	303	AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
	304	}
	305
	306	fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
	307	fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
	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);
	313	fContrH->Write(fContrH->GetName(),TObject::kOverwrite);
	314	fContrHa->Write(fContrHa->GetName(),TObject::kOverwrite);
	315	fmv.Close();
	316	}
	317
	318	//______________________________________________________________________
	319	Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,UInt_t mult){
	320	// Apply selection criteria to events
	321	Bool_t status = kFALSE;
	322	if(!vert)return status;
	323	// Remove vertices reconstructed with vertexerZ
	324	if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
	325	Int_t ncontr = vert->GetNContributors();
	326	AliDebug(1,Form("Number of contributors = %d",ncontr));
	327	Double_t ex = vert->GetXRes();
	328	fMultH->Fill(mult);
	329	fErrXH->Fill(ex*1000.);
	330	fContrH->Fill((Float_t)ncontr);
	331	if(ncontr>fFilterOnContributors && (mult>fLowSPD0 && mult<fHighSPD0) &&
	332	((ex*1000.)<fErrXCut)) {
	333	status = kTRUE;
	334	fMultHa->Fill(mult);
	335	fErrXHa->Fill(ex*1000.);
	336	fContrHa->Fill((Float_t)ncontr);
	337	}
	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++){
	351	fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
	352	fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
	353	fAverPosSum[i]+=currentPos[i];
	354	}
	355	for(Int_t i=0;i<3;i++){
	356	for(Int_t j=i;j<3;j++){
	357	fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
	358	}
	359	}
	360
	361	fTotContributors+=vert->GetNContributors();
	362	fVertexXY->Fill(currentPos[0],currentPos[1]);
	363	fVertexZ->Fill(currentPos[2]);
	364
	365	fNoEventsContr++;
	366	}
	367
	368	//______________________________________________________________________
	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
	373
	374	static Bool_t preliminary = kTRUE;
	375	static Int_t oldnumbevt = 0;
	376	if(!(preliminary \|\| killOutliers))return kTRUE; //ComputeMean is
	377	// executed only once with argument kFALSE
	378	Double_t wpos[3];
	379	for(Int_t i=0;i<3;i++)wpos[i]=fWeighPos[i];
	380
	381	Int_t istart = oldnumbevt;
	382	if(killOutliers)istart = 0;
	383	if(killOutliers && preliminary){
	384	preliminary = kFALSE;
	385	Reset(kTRUE,kFALSE);
	386	}
	387	oldnumbevt = fVertArray.GetEntries();
	388	for(Int_t i=istart;i<fVertArray.GetEntries();i++){
	389	AliESDVertex* vert = NULL;
	390
	391	// second pass
	392	if(killOutliers && fAccEvents.TestBitNumber(i)){
	393	vert=(AliESDVertex*)fVertArray[i];
	394	Double_t dist=(vert->GetXv()-wpos[0])*(vert->GetXv()-wpos[0]);
	395	dist+=(vert->GetYv()-wpos[1])*(vert->GetYv()-wpos[1]);
	396	dist=sqrt(dist)*10.; // distance in mm
	397	fDistH->Fill(dist);
	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	}
	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	}
	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	}
	430	if(killOutliers)ResetArray();
	431	fAverContributors = fTotContributors/nevents;
	432	return kTRUE;
	433	}