[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 "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 "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(Bool_t mode):TObject(),
fDetTypeRec(NULL),
fVertexXY(NULL),
fVertexZ(NULL),
fNoEventsContr(0),
fTotTracklets(0.),
fAverTracklets(0.),
fTotTrackletsSq(0.),
fSigmaOnAverTracks(0.), 
fFilterOnContributors(0),
fFilterOnTracklets(0),
fMode(mode),
fVertexer(NULL)
{
  // 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);
  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");
}

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

  // Instatiate vertexer
  if (!fMode) {
    fVertexer = new AliITSVertexer3DTapan(1000);
  }
  else {
    fVertexer = new AliITSVertexer3D();
    fVertexer->SetDetTypeRec(fDetTypeRec);
    fVertexer->SetComputeMultiplicity(kTRUE);
  }
  return kTRUE;
}

//______________________________________________________________________
AliITSMeanVertexer::~AliITSMeanVertexer() {
  // Destructor
  delete fDetTypeRec;
  delete fVertexXY;
  delete fVertexZ;
  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
  AliITSRecPointContainer::Instance()->PrepareToRead();
  TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
  fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");

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