[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 = 2.0, yLimit = 2.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;
  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
a3a3a28f	1	#include <TFile.h>
37add1d1	2	#include <TH1.h>
37add1d1	3	#include <TH2.h>
a3a3a28f	4	#include "AliGeomManager.h"
a3a3a28f	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"
37add1d1	15	#include "AliITSVertexer3DTapan.h"
a3a3a28f	16	#include "AliESDVertex.h"
	17	#include "AliMeanVertex.h"
	18	#include "AliMultiplicity.h"
	19
a3a3a28f	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(),
37add1d1	42	fDetTypeRec(NULL),
	43	fVertexXY(NULL),
	44	fVertexZ(NULL),
a3a3a28f	45	fNoEventsContr(0),
	46	fTotTracklets(0.),
	47	fAverTracklets(0.),
762133f4	48	fTotTrackletsSq(0.),
a3a3a28f	49	fSigmaOnAverTracks(0.),
a3a3a28f	50	fFilterOnContributors(0),
37add1d1	51	fFilterOnTracklets(0)
a3a3a28f	52	{
	53	// Default Constructor
	54	for(Int_t i=0;i<3;i++){
762133f4	55	fWeighPosSum[i] = 0.;
	56	fWeighSigSum[i] = 0.;
	57	fAverPosSum[i] = 0.;
a3a3a28f	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.;
762133f4	62	for(Int_t j=0; j<3;j++)fAverPosSqSum[i][j] = 0.;
a3a3a28f	63	}
a3a3a28f	64
37add1d1	65	// Histograms initialization
	66	const Float_t xLimit = 2.0, yLimit = 2.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);
a3a3a28f	73	}
	74
	75	//______________________________________________________________________
37add1d1	76	Bool_t AliITSMeanVertexer::Init() {
37add1d1	77	// Initialize filters
a3a3a28f	78	// Initialize geometry
37add1d1	79	// Initialize ITS classes
a3a3a28f	80
37add1d1	81	AliGeomManager::LoadGeometry();
37add1d1	82	if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
a3a3a28f	83
	84	AliITSInitGeometry initgeom;
	85	AliITSgeom *geom = initgeom.CreateAliITSgeom();
37add1d1	86	if (!geom) return kFALSE;
a3a3a28f	87	printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
37add1d1	88
cf0b1cea	89	fDetTypeRec = new AliITSDetTypeRec();
762133f4	90	fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
37add1d1	91	fDetTypeRec->SetITSgeom(geom);
	92	fDetTypeRec->SetDefaults();
	93	fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
	94
a3a3a28f	95	// Initialize filter values to their defaults
	96	SetFilterOnContributors();
	97	SetFilterOnTracklets();
a3a3a28f	98
37add1d1	99	return kTRUE;
a3a3a28f	100	}
	101
	102	//______________________________________________________________________
	103	AliITSMeanVertexer::~AliITSMeanVertexer() {
	104	// Destructor
37add1d1	105	delete fDetTypeRec;
	106	delete fVertexXY;
	107	delete fVertexZ;
a3a3a28f	108	}
	109
	110	//______________________________________________________________________
37add1d1	111	Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader, Bool_t mode){
a3a3a28f	112	// Performs SPD local reconstruction
37add1d1	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;
a3a3a28f	128	}
37add1d1	129	else {
	130	// Run standard vertexer3d
	131	AliITSVertexer3D *vertexer2 = new AliITSVertexer3D();
cf0b1cea	132	vtx = vertexer2->FindVertexForCurrentEvent(clustersTree);
37add1d1	133	AliMultiplicity *mult = vertexer2->GetMultiplicity();
	134	delete vertexer2;
	135	if(Filter(vtx,mult)) vtxOK = kTRUE;
a3a3a28f	136	}
37add1d1	137	delete clustersTree;
	138	if (vtxOK) AddToMean(vtx);
	139	if (vtx) delete vtx;
a3a3a28f	140
37add1d1	141	return vtxOK;
a3a3a28f	142	}
	143
	144	//______________________________________________________________________
37add1d1	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
a3a3a28f	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
37add1d1	160	AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,fTotTracklets,fAverTracklets,fSigmaOnAverTracks);
a3a3a28f	161	mv.SetTitle("Mean Vertex");
37add1d1	162	mv.SetName("MeanVertex");
a3a3a28f	163	AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets()));
a3a3a28f	164	AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks()));
37add1d1	165	// we have to add chi2 here
762133f4	166	AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverTracklets),"MeanVertexPos");
37add1d1	167
	168	mv.Write(mv.GetName(),TObject::kOverwrite);
	169	vtx.Write(vtx.GetName(),TObject::kOverwrite);
a3a3a28f	170	}
	171	else {
	172	AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
	173	}
37add1d1	174
	175	fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
	176	fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
	177	fmv.Close();
a3a3a28f	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	Int_t ncontr = vert->GetNContributors();
	186	Int_t ntracklets = mult->GetNumberOfTracklets();
	187	AliDebug(1,Form("Number of contributors = %d",ncontr));
	188	AliDebug(1,Form("Number of tracklets = %d",ntracklets));
	189	if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE;
762133f4	190	fTotTracklets += ntracklets;
762133f4	191	fTotTrackletsSq += ntracklets*ntracklets;
a3a3a28f	192	return status;
	193	}
	194
	195	//______________________________________________________________________
	196	void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
	197	// update mean vertex
	198	Double_t currentPos[3],currentSigma[3];
	199	vert->GetXYZ(currentPos);
	200	vert->GetSigmaXYZ(currentSigma);
	201	Bool_t goon = kTRUE;
	202	for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
	203	if(!goon)return;
	204	for(Int_t i=0;i<3;i++){
762133f4	205	fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
	206	fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
	207	fAverPosSum[i]+=currentPos[i];
a3a3a28f	208	}
	209	for(Int_t i=0;i<3;i++){
	210	for(Int_t j=i;j<3;j++){
762133f4	211	fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
a3a3a28f	212	}
a3a3a28f	213	}
37add1d1	214
	215	fVertexXY->Fill(currentPos[0],currentPos[1]);
	216	fVertexZ->Fill(currentPos[2]);
	217
a3a3a28f	218	fNoEventsContr++;
	219	}
	220
	221	//______________________________________________________________________
	222	Bool_t AliITSMeanVertexer::ComputeMean(){
	223	// compute mean vertex
	224	if(fNoEventsContr < 2) return kFALSE;
	225	Double_t nevents = fNoEventsContr;
	226	for(Int_t i=0;i<3;i++){
762133f4	227	fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i];
	228	fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
	229	fAverPos[i] = fAverPosSum[i]/nevents;
a3a3a28f	230	}
	231	for(Int_t i=0;i<3;i++){
	232	for(Int_t j=i;j<3;j++){
762133f4	233	fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
a3a3a28f	234	fAverPosSq[i][j] -= nevents/(nevents -1.)fAverPos[i]fAverPos[j];
	235	}
	236	}
762133f4	237
	238	fAverTracklets = fTotTracklets/nevents;
	239	fSigmaOnAverTracks = fTotTrackletsSq/(nevents - 1);
	240	fSigmaOnAverTracks -= nevents/(nevents -1.)fAverTrackletsfAverTracklets;
a3a3a28f	241	fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks);
	242	return kTRUE;
	243	}