[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);
  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();

  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();
    vertexer2->SetDetTypeRec(fDetTypeRec);
    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
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
8d7f6205	66	const Float_t xLimit = 5.0, yLimit = 5.0, zLimit = 50.0;
37add1d1	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);
63845149	71	fVertexXY->SetXTitle("X , cm");
	72	fVertexXY->SetYTitle("Y , cm");
	73	fVertexXY->SetOption("colz");
37add1d1	74	fVertexZ = new TH1F("VertexZ"," Longitudinal Vertex Profile",
37add1d1	75	2*(Int_t)(zLimit/zDelta),-zLimit,zLimit);
63845149	76	fVertexZ->SetXTitle("Z , cm");
a3a3a28f	77	}
	78
	79	//______________________________________________________________________
37add1d1	80	Bool_t AliITSMeanVertexer::Init() {
37add1d1	81	// Initialize filters
a3a3a28f	82	// Initialize geometry
37add1d1	83	// Initialize ITS classes
a3a3a28f	84
37add1d1	85	AliGeomManager::LoadGeometry();
37add1d1	86	if (!AliGeomManager::ApplyAlignObjsFromCDB("ITS")) return kFALSE;
a3a3a28f	87
	88	AliITSInitGeometry initgeom;
	89	AliITSgeom *geom = initgeom.CreateAliITSgeom();
37add1d1	90	if (!geom) return kFALSE;
a3a3a28f	91	printf("Geometry name: %s \n",(initgeom.GetGeometryName()).Data());
37add1d1	92
cf0b1cea	93	fDetTypeRec = new AliITSDetTypeRec();
762133f4	94	fDetTypeRec->SetLoadOnlySPDCalib(kTRUE);
37add1d1	95	fDetTypeRec->SetITSgeom(geom);
	96	fDetTypeRec->SetDefaults();
	97	fDetTypeRec->SetDefaultClusterFindersV2(kTRUE);
	98
a3a3a28f	99	// Initialize filter values to their defaults
	100	SetFilterOnContributors();
	101	SetFilterOnTracklets();
a3a3a28f	102
37add1d1	103	return kTRUE;
a3a3a28f	104	}
	105
	106	//______________________________________________________________________
	107	AliITSMeanVertexer::~AliITSMeanVertexer() {
	108	// Destructor
37add1d1	109	delete fDetTypeRec;
	110	delete fVertexXY;
	111	delete fVertexZ;
a3a3a28f	112	}
	113
	114	//______________________________________________________________________
37add1d1	115	Bool_t AliITSMeanVertexer::Reconstruct(AliRawReader *rawReader, Bool_t mode){
a3a3a28f	116	// Performs SPD local reconstruction
37add1d1	117	// and vertex finding
	118	// returns true in case a vertex is found
	119
	120	// Run SPD cluster finder
	121	TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
	122	fDetTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
	123
	124	Bool_t vtxOK = kFALSE;
	125	AliESDVertex *vtx = NULL;
	126	// Run Tapan's vertex-finder
	127	if (!mode) {
	128	AliITSVertexer3DTapan *vertexer1 = new AliITSVertexer3DTapan(1000);
	129	vtx = vertexer1->FindVertexForCurrentEvent(clustersTree);
	130	delete vertexer1;
	131	if (TMath::Abs(vtx->GetChi2()) < 0.1) vtxOK = kTRUE;
a3a3a28f	132	}
37add1d1	133	else {
	134	// Run standard vertexer3d
	135	AliITSVertexer3D *vertexer2 = new AliITSVertexer3D();
307ebfc0	136	vertexer2->SetDetTypeRec(fDetTypeRec);
cf0b1cea	137	vtx = vertexer2->FindVertexForCurrentEvent(clustersTree);
37add1d1	138	AliMultiplicity *mult = vertexer2->GetMultiplicity();
	139	delete vertexer2;
	140	if(Filter(vtx,mult)) vtxOK = kTRUE;
a3a3a28f	141	}
37add1d1	142	delete clustersTree;
	143	if (vtxOK) AddToMean(vtx);
	144	if (vtx) delete vtx;
a3a3a28f	145
37add1d1	146	return vtxOK;
a3a3a28f	147	}
	148
	149	//______________________________________________________________________
37add1d1	150	void AliITSMeanVertexer::WriteVertices(const char *filename){
	151	// Compute mean vertex and
	152	// store it along with the histograms
	153	// in a file
	154
	155	TFile fmv(filename,"update");
	156
a3a3a28f	157	if(ComputeMean()){
	158	Double_t cov[6];
	159	cov[0] = fAverPosSq[0][0]; // variance x
	160	cov[1] = fAverPosSq[0][1]; // cov xy
	161	cov[2] = fAverPosSq[1][1]; // variance y
	162	cov[3] = fAverPosSq[0][2]; // cov xz
	163	cov[4] = fAverPosSq[1][2]; // cov yz
	164	cov[5] = fAverPosSq[2][2]; // variance z
37add1d1	165	AliMeanVertex mv(fWeighPos,fWeighSig,cov,fNoEventsContr,fTotTracklets,fAverTracklets,fSigmaOnAverTracks);
a3a3a28f	166	mv.SetTitle("Mean Vertex");
37add1d1	167	mv.SetName("MeanVertex");
a3a3a28f	168	AliDebug(1,Form("Contrib av. trk = %10.2f ",mv.GetAverageNumbOfTracklets()));
a3a3a28f	169	AliDebug(1,Form("Sigma %10.4f ",mv.GetSigmaOnAvNumbOfTracks()));
37add1d1	170	// we have to add chi2 here
762133f4	171	AliESDVertex vtx(fWeighPos,cov,0,TMath::Nint(fAverTracklets),"MeanVertexPos");
37add1d1	172
	173	mv.Write(mv.GetName(),TObject::kOverwrite);
	174	vtx.Write(vtx.GetName(),TObject::kOverwrite);
a3a3a28f	175	}
	176	else {
	177	AliError(Form("Evaluation of mean vertex not possible. Number of used events = %d",fNoEventsContr));
	178	}
37add1d1	179
	180	fVertexXY->Write(fVertexXY->GetName(),TObject::kOverwrite);
	181	fVertexZ->Write(fVertexZ->GetName(),TObject::kOverwrite);
	182	fmv.Close();
a3a3a28f	183	}
	184
	185	//______________________________________________________________________
	186	Bool_t AliITSMeanVertexer::Filter(AliESDVertex vert,AliMultiplicity mult){
	187	// Apply selection criteria to events
	188	Bool_t status = kFALSE;
	189	if(!vert \|\| !mult)return status;
8d7f6205	190	// Remove vertices reconstructed with vertexerZ
8d7f6205	191	if(strcmp(vert->GetName(),"SPDVertexZ") == 0) return status;
a3a3a28f	192	Int_t ncontr = vert->GetNContributors();
	193	Int_t ntracklets = mult->GetNumberOfTracklets();
	194	AliDebug(1,Form("Number of contributors = %d",ncontr));
	195	AliDebug(1,Form("Number of tracklets = %d",ntracklets));
	196	if(ncontr>fFilterOnContributors && ntracklets > fFilterOnTracklets) status = kTRUE;
762133f4	197	fTotTracklets += ntracklets;
762133f4	198	fTotTrackletsSq += ntracklets*ntracklets;
a3a3a28f	199	return status;
	200	}
	201
	202	//______________________________________________________________________
	203	void AliITSMeanVertexer::AddToMean(AliESDVertex *vert){
	204	// update mean vertex
	205	Double_t currentPos[3],currentSigma[3];
	206	vert->GetXYZ(currentPos);
	207	vert->GetSigmaXYZ(currentSigma);
	208	Bool_t goon = kTRUE;
	209	for(Int_t i=0;i<3;i++)if(currentSigma[i] == 0.)goon = kFALSE;
	210	if(!goon)return;
	211	for(Int_t i=0;i<3;i++){
762133f4	212	fWeighPosSum[i]+=currentPos[i]/currentSigma[i]/currentSigma[i];
	213	fWeighSigSum[i]+=1./currentSigma[i]/currentSigma[i];
	214	fAverPosSum[i]+=currentPos[i];
a3a3a28f	215	}
	216	for(Int_t i=0;i<3;i++){
	217	for(Int_t j=i;j<3;j++){
762133f4	218	fAverPosSqSum[i][j] += currentPos[i] * currentPos[j];
a3a3a28f	219	}
a3a3a28f	220	}
37add1d1	221
	222	fVertexXY->Fill(currentPos[0],currentPos[1]);
	223	fVertexZ->Fill(currentPos[2]);
	224
a3a3a28f	225	fNoEventsContr++;
	226	}
	227
	228	//______________________________________________________________________
	229	Bool_t AliITSMeanVertexer::ComputeMean(){
	230	// compute mean vertex
	231	if(fNoEventsContr < 2) return kFALSE;
	232	Double_t nevents = fNoEventsContr;
	233	for(Int_t i=0;i<3;i++){
762133f4	234	fWeighPos[i] = fWeighPosSum[i]/fWeighSigSum[i];
	235	fWeighSig[i] = 1./TMath::Sqrt(fWeighSigSum[i]);
	236	fAverPos[i] = fAverPosSum[i]/nevents;
a3a3a28f	237	}
	238	for(Int_t i=0;i<3;i++){
	239	for(Int_t j=i;j<3;j++){
762133f4	240	fAverPosSq[i][j] = fAverPosSqSum[i][j]/(nevents -1.);
a3a3a28f	241	fAverPosSq[i][j] -= nevents/(nevents -1.)fAverPos[i]fAverPos[j];
	242	}
	243	}
762133f4	244
	245	fAverTracklets = fTotTracklets/nevents;
	246	fSigmaOnAverTracks = fTotTrackletsSq/(nevents - 1);
	247	fSigmaOnAverTracks -= nevents/(nevents -1.)fAverTrackletsfAverTracklets;
a3a3a28f	248	fSigmaOnAverTracks = TMath::Sqrt(fSigmaOnAverTracks);
	249	return kTRUE;
	250	}