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

/**************************************************************************
 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

#include <TClonesArray.h>
#include "AliLog.h"
#include "AliESDVertex.h"
#include "AliRunLoader.h"
#include "AliITSLoader.h"
#include "AliITSgeomTGeo.h"
#include "AliITSRecPoint.h"
#include "AliITSVertexerCosmics.h"
#include "AliStrLine.h"

//------------------------------------------------------------------------
// This class implements a method to construct a "fake" primary
// vertex for cosmic events in which the muon crosses one of 5 inner
// ITS layers. A fake primary vertex is needed for the reconstruction,
// with e.g. AliITStrackerSA, of the two tracks produced by the muon 
// in the ITS.
//   We build pairs of clusters on a given layer and define the fake vertex as
// the mid-point of the straight line joining the two clusters.
//   We use the innermost layer that has at least two clusters.
//   We reject the background by requiring at least one cluster on the outer
// layer, closer than fMaxDistOnOuterLayer to the tracklet prolongation.
//   We can reject (potentially pathological) events with the muon track
// tangential to the layer by the requiring the radial position of
// the vertex to be smaller than fMaxVtxRadius.
//   Due to background clusters, more than one vertex per event can 
// be found. We consider the first found.
//   The errors on x,y,z of the vertex are calculated as errors on the mean
// of clusters coordinates. Non-diag elements of vertex cov. mat. are set to 0.
//   The number of contributors set in the AliESDVertex object is the
// number of the layer on which the tracklet was built; if this number is -1, 
// the procedure could not find a vertex position and by default 
// the vertex coordinates are set to (0,0,0) with large errors (100,100,100)
//
// Origin: A.Dainese, andrea.dainese@lnl.infn.it
//-------------------------------------------------------------------------

ClassImp(AliITSVertexerCosmics)

//-------------------------------------------------------------------------
AliITSVertexerCosmics::AliITSVertexerCosmics():AliITSVertexer(),
fMaxDistOnOuterLayer(0),
fMinDist2Vtxs(0)
{
  // Default constructor
  SetFirstLastModules(0,0,79);
  SetFirstLastModules(1,80,239);
  SetFirstLastModules(2,240,323);
  SetFirstLastModules(3,324,499);
  SetFirstLastModules(4,500,1247);
  SetFirstLastModules(5,1248,2197);
  SetMaxVtxRadius(0,3.5);
  SetMaxVtxRadius(1,6.5);
  SetMaxVtxRadius(2,14.5);
  SetMaxVtxRadius(3,23.5);
  SetMaxVtxRadius(4,37.5);
  SetMaxVtxRadius(5,42.5);
  SetMaxDistOnOuterLayer();
  SetMinDist2Vtxs();
}
//--------------------------------------------------------------------------
AliESDVertex* AliITSVertexerCosmics::FindVertexForCurrentEvent(Int_t evnumber) 
{
  // Defines the AliESDVertex for the current event

  fCurrentVertex = 0;
  AliRunLoader *rl =AliRunLoader::GetRunLoader();
  AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
  itsLoader->LoadRecPoints();
  rl->GetEvent(evnumber);

  TTree *rpTree = itsLoader->TreeR();

  TClonesArray *recpoints=new TClonesArray("AliITSRecPoint",10000);
  rpTree->SetBranchAddress("ITSRecPoints",&recpoints);

  Int_t lay,lad,det; 

  // Search for innermost layer with at least two clusters 
  // on two different modules
  Int_t ilayer=0;
  while(ilayer<6) {
    Int_t nHitModules=0;
    for(Int_t imodule=fFirst[ilayer]; imodule<=fLast[ilayer]; imodule++) {
      rpTree->GetEvent(imodule);
      AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
      lay -= 1;  // AliITSgeomTGeo gives layer from 1 to 6, we want 0 to 5
      if(lay!=ilayer) AliFatal("Layer mismatch!");
      if(recpoints->GetEntriesFast()>0) nHitModules++;
    }
    if(nHitModules>=2) break;
    ilayer++;
  }
  printf("Building tracklets on layer %d\n",ilayer);

  const Int_t arrSize = 1000;
  Float_t xclInnLay[arrSize],yclInnLay[arrSize],zclInnLay[arrSize],modclInnLay[arrSize];
  Float_t e2xclInnLay[arrSize],e2yclInnLay[arrSize],e2zclInnLay[arrSize];
  Int_t nclInnLayStored=0;
  Float_t xclOutLay[arrSize],yclOutLay[arrSize],zclOutLay[arrSize],modclOutLay[arrSize];
  Int_t nclOutLayStored=0;
  Int_t nRecPoints,nRecPointsInnLay=0;

  Float_t gc[3],gcov[5];

  Float_t x[arrSize],y[arrSize],z[arrSize],e2x[arrSize],e2y[arrSize],e2z[arrSize];
  Double_t p1[3],p2[3],p3[3];
  Int_t nvtxs;
  Bool_t good,matchtoOutLay;
  Float_t xvtx,yvtx,zvtx,rvtx;

  Double_t pos[3]={0.,0.,0.};
  Double_t err[3]={100.,100.,100.};
  Int_t ncontributors = -1;

  // Collect clusters in the selected layer and the outer one
  for(Int_t imodule=fFirst[ilayer]; imodule<=fLast[ilayer+1]; imodule++) {
    rpTree->GetEvent(imodule);
    AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
    lay -= 1; // AliITSgeomTGeo gives layer from 1 to 6, we want 0 to 5
    nRecPoints=recpoints->GetEntriesFast();
    if(imodule<=fLast[ilayer]) nRecPointsInnLay += nRecPoints;
    //printf("cosmics: module %d clusters %d\n",imodule,nRecPoints);
    for(Int_t irp=0; irp<nRecPoints; irp++) {
      AliITSRecPoint *rp=(AliITSRecPoint*)recpoints->UncheckedAt(irp);
      // Local coordinates of this recpoint
      rp->GetGlobalXYZ(gc);
      if(lay==ilayer) { // store InnLay clusters
	xclInnLay[nclInnLayStored]=gc[0];
	yclInnLay[nclInnLayStored]=gc[1];
	zclInnLay[nclInnLayStored]=gc[2];
	rp->GetGlobalCov(gcov);
	e2xclInnLay[nclInnLayStored]=gcov[0];
	e2yclInnLay[nclInnLayStored]=gcov[3];
	e2zclInnLay[nclInnLayStored]=gcov[5];
	modclInnLay[nclInnLayStored]=imodule;
	nclInnLayStored++;
      }
      if(lay==ilayer+1) { // store OutLay clusters
	xclOutLay[nclOutLayStored]=gc[0];
	yclOutLay[nclOutLayStored]=gc[1];
	zclOutLay[nclOutLayStored]=gc[2];
	modclOutLay[nclOutLayStored]=imodule;
	nclOutLayStored++;
      }
      if(nclInnLayStored>arrSize || nclOutLayStored>arrSize) {
	//AliFatal("More than arrSize clusters per layer");
	AliWarning("Too many clusters per layer");
	delete recpoints;
	itsLoader->UnloadRecPoints();
	fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
	fCurrentVertex->SetTitle("cosmics fake vertex (failed)");
	fCurrentVertex->SetNContributors(ncontributors);
	return fCurrentVertex;
      }
    }// end clusters in a module
  }// end modules

  // build fake vertices
  nvtxs=0;
  // InnLay - first cluster
  for(Int_t i1InnLay=0; i1InnLay<nclInnLayStored; i1InnLay++) { 
    p1[0]=xclInnLay[i1InnLay]; 
    p1[1]=yclInnLay[i1InnLay]; 
    p1[2]=zclInnLay[i1InnLay];
    // InnLay - second cluster
    for(Int_t i2InnLay=i1InnLay+1; i2InnLay<nclInnLayStored; i2InnLay++) { 
      if(modclInnLay[i1InnLay]==modclInnLay[i2InnLay]) continue;
      p2[0]=xclInnLay[i2InnLay]; 
      p2[1]=yclInnLay[i2InnLay]; 
      p2[2]=zclInnLay[i2InnLay];
      // look for point on OutLay
      AliStrLine InnLayline(p1,p2,kTRUE);
      matchtoOutLay = kFALSE;
      for(Int_t iOutLay=0; iOutLay<nclOutLayStored; iOutLay++) {
	p3[0]=xclOutLay[iOutLay]; 
	p3[1]=yclOutLay[iOutLay]; 
	p3[2]=zclOutLay[iOutLay];
	//printf(" %f\n",InnLayline.GetDistFromPoint(p3));
	if(InnLayline.GetDistFromPoint(p3)<fMaxDistOnOuterLayer) 
	  { matchtoOutLay = kTRUE; break; }
      }
      if(!matchtoOutLay) continue;
      xvtx = 0.5*(xclInnLay[i1InnLay]+xclInnLay[i2InnLay]);
      yvtx = 0.5*(yclInnLay[i1InnLay]+yclInnLay[i2InnLay]);
      zvtx = 0.5*(zclInnLay[i1InnLay]+zclInnLay[i2InnLay]);
      rvtx = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx);
      if(rvtx>fMaxVtxRadius[ilayer]) continue;
      good = kTRUE;
      for(Int_t iv=0; iv<nvtxs; iv++) {
	if(TMath::Sqrt((xvtx- x[iv])*(xvtx- x[iv])+
		       (yvtx- y[iv])*(yvtx- y[iv])+
		       (zvtx- z[iv])*(zvtx- z[iv])) < fMinDist2Vtxs) 
	  good = kFALSE;
      }
      if(good) {
	x[nvtxs]=xvtx;
	y[nvtxs]=yvtx;
	z[nvtxs]=zvtx;
	e2x[nvtxs]=0.25*(e2xclInnLay[i1InnLay]+e2xclInnLay[i2InnLay]);
	e2y[nvtxs]=0.25*(e2yclInnLay[i1InnLay]+e2yclInnLay[i2InnLay]);
	e2z[nvtxs]=0.25*(e2zclInnLay[i1InnLay]+e2zclInnLay[i2InnLay]);
	nvtxs++;
      }
    } // InnLay - second cluster
  } // InnLay - first cluster

  if(nvtxs) { 
    ncontributors = ilayer;
    pos[0]=x[0]; 
    pos[1]=y[0]; 
    pos[2]=z[0];
    err[0]=TMath::Sqrt(e2x[0]); 
    err[1]=TMath::Sqrt(e2y[0]); 
    err[2]=TMath::Sqrt(e2z[0]);
  }

  fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
  fCurrentVertex->SetTitle("cosmics fake vertex");
  fCurrentVertex->SetNContributors(ncontributors);
  //fCurrentVertex->Print();

  delete recpoints;
  itsLoader->UnloadRecPoints();

  return fCurrentVertex;
}  
//-------------------------------------------------------------------------
void AliITSVertexerCosmics::FindVertices()
{
  // computes the vertices of the events in the range FirstEvent - LastEvent
  AliRunLoader *rl = AliRunLoader::GetRunLoader();
  AliITSLoader* itsLoader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
  itsLoader->ReloadRecPoints();
  for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
    //  printf("Processing event %d\n",i);
    rl->GetEvent(i);
    FindVertexForCurrentEvent(i);
    if(fCurrentVertex){
      WriteCurrentVertex();
    }
  }
}
//-------------------------------------------------------------------------
void AliITSVertexerCosmics::PrintStatus() const 
{
  // Print current status
  printf("=======================================================\n");
  printf(" fMaxDistOnOuterLayer: %f\n",fMaxDistOnOuterLayer);
  printf(" fMaxVtxRadius[0]:  %f\n",fMaxVtxRadius[0]);
  printf(" fMinDist2Vtxs:  %f\n",fMinDist2Vtxs);
  printf("=======================================================\n");
}
//-------------------------------------------------------------------------
Commit	Line	Data
3acc14d5	1	/**************************************************************************
	2	* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	#include <TClonesArray.h>
	17	#include "AliLog.h"
	18	#include "AliESDVertex.h"
	19	#include "AliRunLoader.h"
	20	#include "AliITSLoader.h"
0333e459	21	#include "AliITSgeomTGeo.h"
3acc14d5	22	#include "AliITSRecPoint.h"
	23	#include "AliITSVertexerCosmics.h"
	24	#include "AliStrLine.h"
	25
	26	//------------------------------------------------------------------------
	27	// This class implements a method to construct a "fake" primary
b8ed1a92	28	// vertex for cosmic events in which the muon crosses one of 5 inner
b8ed1a92	29	// ITS layers. A fake primary vertex is needed for the reconstruction,
3acc14d5	30	// with e.g. AliITStrackerSA, of the two tracks produced by the muon
3acc14d5	31	// in the ITS.
b8ed1a92	32	// We build pairs of clusters on a given layer and define the fake vertex as
3acc14d5	33	// the mid-point of the straight line joining the two clusters.
b8ed1a92	34	// We use the innermost layer that has at least two clusters.
	35	// We reject the background by requiring at least one cluster on the outer
	36	// layer, closer than fMaxDistOnOuterLayer to the tracklet prolongation.
3acc14d5	37	// We can reject (potentially pathological) events with the muon track
b8ed1a92	38	// tangential to the layer by the requiring the radial position of
3acc14d5	39	// the vertex to be smaller than fMaxVtxRadius.
	40	// Due to background clusters, more than one vertex per event can
	41	// be found. We consider the first found.
b8ed1a92	42	// The errors on x,y,z of the vertex are calculated as errors on the mean
b8ed1a92	43	// of clusters coordinates. Non-diag elements of vertex cov. mat. are set to 0.
3acc14d5	44	// The number of contributors set in the AliESDVertex object is the
b8ed1a92	45	// number of the layer on which the tracklet was built; if this number is -1,
3acc14d5	46	// the procedure could not find a vertex position and by default
3acc14d5	47	// the vertex coordinates are set to (0,0,0) with large errors (100,100,100)
3acc14d5	48	//
	49	// Origin: A.Dainese, andrea.dainese@lnl.infn.it
	50	//-------------------------------------------------------------------------
	51
	52	ClassImp(AliITSVertexerCosmics)
	53
	54	//-------------------------------------------------------------------------
	55	AliITSVertexerCosmics::AliITSVertexerCosmics():AliITSVertexer(),
b8ed1a92	56	fMaxDistOnOuterLayer(0),
3acc14d5	57	fMinDist2Vtxs(0)
	58	{
	59	// Default constructor
b8ed1a92	60	SetFirstLastModules(0,0,79);
	61	SetFirstLastModules(1,80,239);
	62	SetFirstLastModules(2,240,323);
	63	SetFirstLastModules(3,324,499);
	64	SetFirstLastModules(4,500,1247);
	65	SetFirstLastModules(5,1248,2197);
	66	SetMaxVtxRadius(0,3.5);
	67	SetMaxVtxRadius(1,6.5);
	68	SetMaxVtxRadius(2,14.5);
	69	SetMaxVtxRadius(3,23.5);
	70	SetMaxVtxRadius(4,37.5);
	71	SetMaxVtxRadius(5,42.5);
	72	SetMaxDistOnOuterLayer();
3acc14d5	73	SetMinDist2Vtxs();
	74	}
	75	//--------------------------------------------------------------------------
	76	AliESDVertex* AliITSVertexerCosmics::FindVertexForCurrentEvent(Int_t evnumber)
	77	{
	78	// Defines the AliESDVertex for the current event
	79
	80	fCurrentVertex = 0;
	81	AliRunLoader *rl =AliRunLoader::GetRunLoader();
	82	AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
3acc14d5	83	itsLoader->LoadRecPoints();
	84	rl->GetEvent(evnumber);
	85
	86	TTree *rpTree = itsLoader->TreeR();
	87
	88	TClonesArray *recpoints=new TClonesArray("AliITSRecPoint",10000);
	89	rpTree->SetBranchAddress("ITSRecPoints",&recpoints);
	90
3acc14d5	91	Int_t lay,lad,det;
3acc14d5	92
b8ed1a92	93	// Search for innermost layer with at least two clusters
	94	// on two different modules
	95	Int_t ilayer=0;
	96	while(ilayer<6) {
	97	Int_t nHitModules=0;
	98	for(Int_t imodule=fFirst[ilayer]; imodule<=fLast[ilayer]; imodule++) {
	99	rpTree->GetEvent(imodule);
	100	AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
	101	lay -= 1; // AliITSgeomTGeo gives layer from 1 to 6, we want 0 to 5
	102	if(lay!=ilayer) AliFatal("Layer mismatch!");
	103	if(recpoints->GetEntriesFast()>0) nHitModules++;
	104	}
	105	if(nHitModules>=2) break;
	106	ilayer++;
	107	}
	108	printf("Building tracklets on layer %d\n",ilayer);
	109
a1fca0e6	110	const Int_t arrSize = 1000;
	111	Float_t xclInnLay[arrSize],yclInnLay[arrSize],zclInnLay[arrSize],modclInnLay[arrSize];
	112	Float_t e2xclInnLay[arrSize],e2yclInnLay[arrSize],e2zclInnLay[arrSize];
b8ed1a92	113	Int_t nclInnLayStored=0;
a1fca0e6	114	Float_t xclOutLay[arrSize],yclOutLay[arrSize],zclOutLay[arrSize],modclOutLay[arrSize];
b8ed1a92	115	Int_t nclOutLayStored=0;
	116	Int_t nRecPoints,nRecPointsInnLay=0;
	117
	118	Float_t gc[3],gcov[5];
	119
a1fca0e6	120	Float_t x[arrSize],y[arrSize],z[arrSize],e2x[arrSize],e2y[arrSize],e2z[arrSize];
b8ed1a92	121	Double_t p1[3],p2[3],p3[3];
3acc14d5	122	Int_t nvtxs;
b8ed1a92	123	Bool_t good,matchtoOutLay;
b8ed1a92	124	Float_t xvtx,yvtx,zvtx,rvtx;
3acc14d5	125
b107ee21	126	Double_t pos[3]={0.,0.,0.};
	127	Double_t err[3]={100.,100.,100.};
	128	Int_t ncontributors = -1;
	129
b8ed1a92	130	// Collect clusters in the selected layer and the outer one
b8ed1a92	131	for(Int_t imodule=fFirst[ilayer]; imodule<=fLast[ilayer+1]; imodule++) {
3acc14d5	132	rpTree->GetEvent(imodule);
0333e459	133	AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
b8ed1a92	134	lay -= 1; // AliITSgeomTGeo gives layer from 1 to 6, we want 0 to 5
	135	nRecPoints=recpoints->GetEntriesFast();
	136	if(imodule<=fLast[ilayer]) nRecPointsInnLay += nRecPoints;
	137	//printf("cosmics: module %d clusters %d\n",imodule,nRecPoints);
	138	for(Int_t irp=0; irp<nRecPoints; irp++) {
3acc14d5	139	AliITSRecPoint rp=(AliITSRecPoint)recpoints->UncheckedAt(irp);
3acc14d5	140	// Local coordinates of this recpoint
b8ed1a92	141	rp->GetGlobalXYZ(gc);
	142	if(lay==ilayer) { // store InnLay clusters
	143	xclInnLay[nclInnLayStored]=gc[0];
	144	yclInnLay[nclInnLayStored]=gc[1];
	145	zclInnLay[nclInnLayStored]=gc[2];
	146	rp->GetGlobalCov(gcov);
	147	e2xclInnLay[nclInnLayStored]=gcov[0];
	148	e2yclInnLay[nclInnLayStored]=gcov[3];
	149	e2zclInnLay[nclInnLayStored]=gcov[5];
	150	modclInnLay[nclInnLayStored]=imodule;
	151	nclInnLayStored++;
3acc14d5	152	}
b8ed1a92	153	if(lay==ilayer+1) { // store OutLay clusters
	154	xclOutLay[nclOutLayStored]=gc[0];
	155	yclOutLay[nclOutLayStored]=gc[1];
	156	zclOutLay[nclOutLayStored]=gc[2];
	157	modclOutLay[nclOutLayStored]=imodule;
	158	nclOutLayStored++;
3acc14d5	159	}
b107ee21	160	if(nclInnLayStored>arrSize \|\| nclOutLayStored>arrSize) {
	161	//AliFatal("More than arrSize clusters per layer");
	162	AliWarning("Too many clusters per layer");
	163	delete recpoints;
	164	itsLoader->UnloadRecPoints();
	165	fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
	166	fCurrentVertex->SetTitle("cosmics fake vertex (failed)");
	167	fCurrentVertex->SetNContributors(ncontributors);
	168	return fCurrentVertex;
	169	}
3acc14d5	170	}// end clusters in a module
b8ed1a92	171	}// end modules
3acc14d5	172
	173	// build fake vertices
	174	nvtxs=0;
b8ed1a92	175	// InnLay - first cluster
	176	for(Int_t i1InnLay=0; i1InnLay<nclInnLayStored; i1InnLay++) {
	177	p1[0]=xclInnLay[i1InnLay];
	178	p1[1]=yclInnLay[i1InnLay];
	179	p1[2]=zclInnLay[i1InnLay];
	180	// InnLay - second cluster
	181	for(Int_t i2InnLay=i1InnLay+1; i2InnLay<nclInnLayStored; i2InnLay++) {
	182	if(modclInnLay[i1InnLay]==modclInnLay[i2InnLay]) continue;
	183	p2[0]=xclInnLay[i2InnLay];
	184	p2[1]=yclInnLay[i2InnLay];
	185	p2[2]=zclInnLay[i2InnLay];
	186	// look for point on OutLay
	187	AliStrLine InnLayline(p1,p2,kTRUE);
	188	matchtoOutLay = kFALSE;
	189	for(Int_t iOutLay=0; iOutLay<nclOutLayStored; iOutLay++) {
	190	p3[0]=xclOutLay[iOutLay];
	191	p3[1]=yclOutLay[iOutLay];
	192	p3[2]=zclOutLay[iOutLay];
	193	//printf(" %f\n",InnLayline.GetDistFromPoint(p3));
	194	if(InnLayline.GetDistFromPoint(p3)<fMaxDistOnOuterLayer)
	195	{ matchtoOutLay = kTRUE; break; }
3acc14d5	196	}
b8ed1a92	197	if(!matchtoOutLay) continue;
	198	xvtx = 0.5*(xclInnLay[i1InnLay]+xclInnLay[i2InnLay]);
	199	yvtx = 0.5*(yclInnLay[i1InnLay]+yclInnLay[i2InnLay]);
	200	zvtx = 0.5*(zclInnLay[i1InnLay]+zclInnLay[i2InnLay]);
3acc14d5	201	rvtx = TMath::Sqrt(xvtxxvtx+yvtxyvtx);
b8ed1a92	202	if(rvtx>fMaxVtxRadius[ilayer]) continue;
3acc14d5	203	good = kTRUE;
	204	for(Int_t iv=0; iv<nvtxs; iv++) {
	205	if(TMath::Sqrt((xvtx- x[iv])*(xvtx- x[iv])+
	206	(yvtx- y[iv])*(yvtx- y[iv])+
	207	(zvtx- z[iv])*(zvtx- z[iv])) < fMinDist2Vtxs)
	208	good = kFALSE;
	209	}
	210	if(good) {
	211	x[nvtxs]=xvtx;
	212	y[nvtxs]=yvtx;
	213	z[nvtxs]=zvtx;
b8ed1a92	214	e2x[nvtxs]=0.25*(e2xclInnLay[i1InnLay]+e2xclInnLay[i2InnLay]);
	215	e2y[nvtxs]=0.25*(e2yclInnLay[i1InnLay]+e2yclInnLay[i2InnLay]);
	216	e2z[nvtxs]=0.25*(e2zclInnLay[i1InnLay]+e2zclInnLay[i2InnLay]);
3acc14d5	217	nvtxs++;
3acc14d5	218	}
b8ed1a92	219	} // InnLay - second cluster
b8ed1a92	220	} // InnLay - first cluster
3acc14d5	221
3acc14d5	222	if(nvtxs) {
b107ee21	223	ncontributors = ilayer;
3acc14d5	224	pos[0]=x[0];
	225	pos[1]=y[0];
	226	pos[2]=z[0];
b8ed1a92	227	err[0]=TMath::Sqrt(e2x[0]);
	228	err[1]=TMath::Sqrt(e2y[0]);
	229	err[2]=TMath::Sqrt(e2z[0]);
3acc14d5	230	}
b107ee21	231
3acc14d5	232	fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
3acc14d5	233	fCurrentVertex->SetTitle("cosmics fake vertex");
b107ee21	234	fCurrentVertex->SetNContributors(ncontributors);
b107ee21	235	//fCurrentVertex->Print();
3acc14d5	236
	237	delete recpoints;
	238	itsLoader->UnloadRecPoints();
	239
	240	return fCurrentVertex;
	241	}
	242	//-------------------------------------------------------------------------
	243	void AliITSVertexerCosmics::FindVertices()
	244	{
	245	// computes the vertices of the events in the range FirstEvent - LastEvent
	246	AliRunLoader *rl = AliRunLoader::GetRunLoader();
	247	AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader");
	248	itsLoader->ReloadRecPoints();
	249	for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
	250	// printf("Processing event %d\n",i);
	251	rl->GetEvent(i);
	252	FindVertexForCurrentEvent(i);
	253	if(fCurrentVertex){
	254	WriteCurrentVertex();
	255	}
	256	}
	257	}
	258	//-------------------------------------------------------------------------
	259	void AliITSVertexerCosmics::PrintStatus() const
	260	{
	261	// Print current status
	262	printf("=======================================================\n");
b8ed1a92	263	printf(" fMaxDistOnOuterLayer: %f\n",fMaxDistOnOuterLayer);
b8ed1a92	264	printf(" fMaxVtxRadius[0]: %f\n",fMaxVtxRadius[0]);
3acc14d5	265	printf(" fMinDist2Vtxs: %f\n",fMinDist2Vtxs);
3acc14d5	266	printf("=======================================================\n");
	267	}
	268	//-------------------------------------------------------------------------