[u/mrichter/AliRoot.git] / ITS / AliITSMeanVertexer.cxx

#include <TFile.h>
#include <TH1.h>
#include <TH2.h>
#include "AliGeomManager.h"
#include "AliITSDetTypeRec.h"
#include "AliITSInitGeometry.h"
#include "AliITSMeanVertexer.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 "AliMultiplicity.h"

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:
// AliITSMeanVertexer mv("RawDataFileName");
// mv.SetGeometryFileName("FileWithGeometry.root"); // def. geometry.root
// ....  optional setters ....
// mv.Reconstruct();  // SPD local reconstruction
// mv.DoVertices(); 
// Resulting AliMeanVertex object is stored in a file named fMVFileName
///////////////////////////////////////////////////////////////////////

/* $Id$ */

//______________________________________________________________________
AliITSMeanVertexer::AliITSMeanVertexer():TObject(),
fDetTypeRec(NULL),
fVertexXY(NULL),
fVertexZ(NULL),
fNoEventsContr(0),
fTotTracklets(0.),
fAverTracklets(0.),
fTotTrackletsSq(0.),
fSigmaOnAverTracks(0.), 
fFilterOnContributors(0),
fFilterOnTracklets(0)
{
  // Default Constructor
  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.;
  }

  // Histograms initialization
  const Float_t xLimit = 5.0, yLimit = 5.0, zLimit = 50.0;
  const Float_t xDelta = 0.02, yDelta = 0.02, 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);
  fVertexZ  = new TH1F("VertexZ"," Longitudinal Vertex Profile",
		       2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
}

//______________________________________________________________________
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();
  SetFilterOnTracklets();

  return kTRUE;
}

//______________________________________________________________________
AliITSMeanVertexer::~AliITSMeanVertexer() {
  // Destructor
  delete fDetTypeRec;
  delete fVertexXY;
  delete fVertexZ;
}

//______________________________________________________________________
Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader, Bool_t mode){
  // Performs SPD local reconstruction
  // and vertex finding
  // returns true in case a vertex is found

  // Run SPD cluster finder
  TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
  fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");

  Bool_t vtxOK = kFALSE;
  AliESDVertex *vtx = NULL;
  // Run Tapan's vertex-finder
  if (!mode) {
    AliITSVertexer3DTapan *vertexer1 = new AliITSVertexer3DTapan(1000);
    vtx = vertexer1->FindVertexForCurrentEvent(clustersTree);
    delete vertexer1;
    if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
  }
  else {
  // Run standard vertexer3d
    AliITSVertexer3D *vertexer2 = new AliITSVertexer3D();
    vtx = vertexer2->FindVertexForCurrentEvent(clustersTree);
    AliMultiplicity *mult = vertexer2->GetMultiplicity();
    delete vertexer2;
    if(Filter(vtx,mult)) vtxOK = kTRUE;
  }
  delete clustersTree;
  if (vtxOK) AddToMean(vtx);
  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");

  if(ComputeMean()){
    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,fTotTracklets,fAverTracklets,fSigmaOnAverTracks);
    mv.SetTitle("Mean Vertex");
    mv.SetName("MeanVertex");
    AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets()));
    AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks()));
    // we have to add chi2 here
    AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverTracklets),"MeanVertexPos");

    mv.Write(mv.GetName(),TObject::kOverwrite);
    vtx.Write(vtx.GetName(),TObject::kOverwrite);
  }
  else {
    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);
  fmv.Close();
}

//______________________________________________________________________
Bool_t AliITSMeanVertexer::Filter(AliESDVertex *vert,AliMultiplicity *mult){
  // Apply selection criteria to events
  Bool_t status = kFALSE;
  if(!vert || !mult)return status;
  // Remove vertices reconstructed with vertexerZ
  if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
  Int_t ncontr = vert->GetNContributors();
  Int_t ntracklets = mult->GetNumberOfTracklets();
  AliDebug(1,Form("Number of contributors = %d",ncontr));
  AliDebug(1,Form("Number of tracklets = %d",ntracklets));
  if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE;
  fTotTracklets += ntracklets;
  fTotTrackletsSq += ntracklets*ntracklets;
  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];
    }
  }

  fVertexXY->Fill(currentPos[0],currentPos[1]);
  fVertexZ->Fill(currentPos[2]);

  fNoEventsContr++;
}

//______________________________________________________________________
Bool_t AliITSMeanVertexer::ComputeMean(){
  // compute mean vertex
  if(fNoEventsContr < 2) return kFALSE;
  Double_t nevents = fNoEventsContr;
  for(Int_t i=0;i<3;i++){
    fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i]; 
    fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
    fAverPos[i] = fAverPosSum[i]/nevents;
  } 
  for(Int_t i=0;i<3;i++){
    for(Int_t j=i;j<3;j++){
      fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
      fAverPosSq[i][j] -= nevents/(nevents -1.)*fAverPos[i]*fAverPos[j];
    }
  } 

  fAverTracklets = fTotTracklets/nevents;
  fSigmaOnAverTracks = fTotTrackletsSq/(nevents - 1);
  fSigmaOnAverTracks -= nevents/(nevents -1.)*fAverTracklets*fAverTracklets;
  fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks);
  return kTRUE;
}
Commit	Line	Data
	1	#include <TFile.h>
	2	#include <TH1.h>
	3	#include <TH2.h>
	4	#include "AliGeomManager.h"
	5	#include "AliITSDetTypeRec.h"
	6	#include "AliITSInitGeometry.h"
	7	#include "AliITSMeanVertexer.h"
	8	#include "AliITSLoader.h"
	9	#include "AliLog.h"
	10	#include "AliRawReader.h"
	11	#include "AliRawReaderDate.h"
	12	#include "AliRawReaderRoot.h"
	13	#include "AliRunLoader.h"
	14	#include "AliITSVertexer3D.h"
	15	#include "AliITSVertexer3DTapan.h"
	16	#include "AliESDVertex.h"
	17	#include "AliMeanVertex.h"
	18	#include "AliMultiplicity.h"
	19
	20	ClassImp(AliITSMeanVertexer)
	21
	22	///////////////////////////////////////////////////////////////////////
	23	// //
	24	// Class to compute vertex position using SPD local reconstruction //
	25	// An average vertex position using all the events //
	26	// is built and saved //
	27	// Individual vertices can be optionally stored //
	28	// Origin: M.Masera (masera@to.infn.it) //
	29	// Usage:
	30	// AliITSMeanVertexer mv("RawDataFileName");
	31	// mv.SetGeometryFileName("FileWithGeometry.root"); // def. geometry.root
	32	// .... optional setters ....
	33	// mv.Reconstruct(); // SPD local reconstruction
	34	// mv.DoVertices();
	35	// Resulting AliMeanVertex object is stored in a file named fMVFileName
	36	///////////////////////////////////////////////////////////////////////
	37
	38	/* $Id$ */
	39
	40	//______________________________________________________________________
	41	AliITSMeanVertexer::AliITSMeanVertexer():TObject(),
	42	fDetTypeRec(NULL),
	43	fVertexXY(NULL),
	44	fVertexZ(NULL),
	45	fNoEventsContr(0),
	46	fTotTracklets(0.),
	47	fAverTracklets(0.),
	48	fTotTrackletsSq(0.),
	49	fSigmaOnAverTracks(0.),
	50	fFilterOnContributors(0),
	51	fFilterOnTracklets(0)
	52	{
	53	// Default Constructor
	54	for(Int_t i=0;i<3;i++){
	55	fWeighPosSum[i] = 0.;
	56	fWeighSigSum[i] = 0.;
	57	fAverPosSum[i] = 0.;
	58	fWeighPos[i] = 0.;
	59	fWeighSig[i] = 0.;
	60	fAverPos[i] = 0.;
	61	for(Int_t j=0; j<3;j++)fAverPosSq[i][j] = 0.;
	62	for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
	63	}
	64
	65	// Histograms initialization
	66	const Float_t xLimit = 5.0, yLimit = 5.0, zLimit = 50.0;
	67	const Float_t xDelta = 0.02, yDelta = 0.02, zDelta = 0.2;
	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);
	71	fVertexZ = new TH1F("VertexZ"," Longitudinal Vertex Profile",
	72	2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
	73	}
	74
	75	//______________________________________________________________________
	76	Bool_t AliITSMeanVertexer::Init() {
	77	// Initialize filters
	78	// Initialize geometry
	79	// Initialize ITS classes
	80
	81	AliGeomManager::LoadGeometry();
	82	if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
	83
	84	AliITSInitGeometry initgeom;
	85	AliITSgeom *geom = initgeom.CreateAliITSgeom();
	86	if (!geom) return kFALSE;
	87	printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
	88
	89	fDetTypeRec = new AliITSDetTypeRec();
	90	fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
	91	fDetTypeRec->SetITSgeom(geom);
	92	fDetTypeRec->SetDefaults();
	93	fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
	94
	95	// Initialize filter values to their defaults
	96	SetFilterOnContributors();
	97	SetFilterOnTracklets();
	98
	99	return kTRUE;
	100	}
	101
	102	//______________________________________________________________________
	103	AliITSMeanVertexer::~AliITSMeanVertexer() {
	104	// Destructor
	105	delete fDetTypeRec;
	106	delete fVertexXY;
	107	delete fVertexZ;
	108	}
	109
	110	//______________________________________________________________________
	111	Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader, Bool_t mode){
	112	// Performs SPD local reconstruction
	113	// and vertex finding
	114	// returns true in case a vertex is found
	115
	116	// Run SPD cluster finder
	117	TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
	118	fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
	119
	120	Bool_t vtxOK = kFALSE;
	121	AliESDVertex *vtx = NULL;
	122	// Run Tapan's vertex-finder
	123	if (!mode) {
	124	AliITSVertexer3DTapan *vertexer1 = new AliITSVertexer3DTapan(1000);
	125	vtx = vertexer1->FindVertexForCurrentEvent(clustersTree);
	126	delete vertexer1;
	127	if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
	128	}
	129	else {
	130	// Run standard vertexer3d
	131	AliITSVertexer3D *vertexer2 = new AliITSVertexer3D();
	132	vtx = vertexer2->FindVertexForCurrentEvent(clustersTree);
	133	AliMultiplicity *mult = vertexer2->GetMultiplicity();
	134	delete vertexer2;
	135	if(Filter(vtx,mult)) vtxOK = kTRUE;
	136	}
	137	delete clustersTree;
	138	if (vtxOK) AddToMean(vtx);
	139	if (vtx) delete vtx;
	140
	141	return vtxOK;
	142	}
	143
	144	//______________________________________________________________________
	145	void AliITSMeanVertexer::WriteVertices(const char *filename){
	146	// Compute mean vertex and
	147	// store it along with the histograms
	148	// in a file
	149
	150	TFile fmv(filename,"update");
	151
	152	if(ComputeMean()){
	153	Double_t cov[6];
	154	cov[0] = fAverPosSq[0][0]; // variance x
	155	cov[1] = fAverPosSq[0][1]; // cov xy
	156	cov[2] = fAverPosSq[1][1]; // variance y
	157	cov[3] = fAverPosSq[0][2]; // cov xz
	158	cov[4] = fAverPosSq[1][2]; // cov yz
	159	cov[5] = fAverPosSq[2][2]; // variance z
	160	AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,fTotTracklets,fAverTracklets,fSigmaOnAverTracks);
	161	mv.SetTitle("Mean Vertex");
	162	mv.SetName("MeanVertex");
	163	AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets()));
	164	AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks()));
	165	// we have to add chi2 here
	166	AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverTracklets),"MeanVertexPos");
	167
	168	mv.Write(mv.GetName(),TObject::kOverwrite);
	169	vtx.Write(vtx.GetName(),TObject::kOverwrite);
	170	}
	171	else {
	172	AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
	173	}
	174
	175	fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
	176	fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
	177	fmv.Close();
	178	}
	179
	180	//______________________________________________________________________
	181	Bool_t AliITSMeanVertexer::Filter(AliESDVertex vert,AliMultiplicity mult){
	182	// Apply selection criteria to events
	183	Bool_t status = kFALSE;
	184	if(!vert \|\| !mult)return status;
	185	// Remove vertices reconstructed with vertexerZ
	186	if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
	187	Int_t ncontr = vert->GetNContributors();
	188	Int_t ntracklets = mult->GetNumberOfTracklets();
	189	AliDebug(1,Form("Number of contributors = %d",ncontr));
	190	AliDebug(1,Form("Number of tracklets = %d",ntracklets));
	191	if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE;
	192	fTotTracklets += ntracklets;
	193	fTotTrackletsSq += ntracklets*ntracklets;
	194	return status;
	195	}
	196
	197	//______________________________________________________________________
	198	void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
	199	// update mean vertex
	200	Double_t currentPos[3],currentSigma[3];
	201	vert->GetXYZ(currentPos);
	202	vert->GetSigmaXYZ(currentSigma);
	203	Bool_t goon = kTRUE;
	204	for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
	205	if(!goon)return;
	206	for(Int_t i=0;i<3;i++){
	207	fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
	208	fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
	209	fAverPosSum[i]+=currentPos[i];
	210	}
	211	for(Int_t i=0;i<3;i++){
	212	for(Int_t j=i;j<3;j++){
	213	fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
	214	}
	215	}
	216
	217	fVertexXY->Fill(currentPos[0],currentPos[1]);
	218	fVertexZ->Fill(currentPos[2]);
	219
	220	fNoEventsContr++;
	221	}
	222
	223	//______________________________________________________________________
	224	Bool_t AliITSMeanVertexer::ComputeMean(){
	225	// compute mean vertex
	226	if(fNoEventsContr < 2) return kFALSE;
	227	Double_t nevents = fNoEventsContr;
	228	for(Int_t i=0;i<3;i++){
	229	fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i];
	230	fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
	231	fAverPos[i] = fAverPosSum[i]/nevents;
	232	}
	233	for(Int_t i=0;i<3;i++){
	234	for(Int_t j=i;j<3;j++){
	235	fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
	236	fAverPosSq[i][j] -= nevents/(nevents -1.)fAverPos[i]fAverPos[j];
	237	}
	238	}
	239
	240	fAverTracklets = fTotTracklets/nevents;
	241	fSigmaOnAverTracks = fTotTrackletsSq/(nevents - 1);
	242	fSigmaOnAverTracks -= nevents/(nevents -1.)fAverTrackletsfAverTracklets;
	243	fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks);
	244	return kTRUE;
	245	}