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

#include <AliITSVertexerPPZ.h>
#include <Riostream.h>
#include <TArrayF.h>
#include <TDirectory.h>
#include <TH1F.h>
#include <TMath.h>
#include <TObjArray.h>
#include <TTree.h>
#include <TClonesArray.h>
#include "AliRun.h"
#include "AliITS.h"
#include "AliITSgeom.h"
#include "AliITSRecPoint.h"
/////////////////////////////////////////////////////////////////////////
//                                                                     //
// This class is intended to compute the Z coordinate of the           //
// primary vertex for p-p interactions, or in general when the         //
// number of secondaries is too slow to use the class                  //
// AliITSVertexerIons, which in turn is optimized for A-A collsions    //
// Origin: masera@to.infn.it      9/12/2002                            //
//                                                                     //
/////////////////////////////////////////////////////////////////////////
ClassImp(AliITSVertexerPPZ)


//______________________________________________________________________
AliITSVertexerPPZ::AliITSVertexerPPZ():AliITSVertexer() {
  // Default Constructor

  SetDiffPhiMax(0);
  fX0 = 0.;
  fY0 = 0.;
  SetFirstLayerModules(0);
  SetSecondLayerModules(0);
  fZFound = 0;
  fZsig = 0.;
  fITS = 0;
}

AliITSVertexerPPZ::AliITSVertexerPPZ(TFile *infile, TFile *outfile, Float_t x0, Float_t y0):AliITSVertexer(infile,outfile) {
  // Standard constructor
  if(!fInFile)Fatal("AliITSVertexerPPZ","No inputfile has been defined!");
  SetDiffPhiMax();
  fX0 = x0;
  fY0 = y0;
  SetFirstLayerModules();
  SetSecondLayerModules();
  fZFound = 0;
  fZsig = 0.;
  fITS = 0;
}

//______________________________________________________________________
AliITSVertexerPPZ::~AliITSVertexerPPZ() {
  // Default Destructor
  fITS = 0;
}

//________________________________________________________
void  AliITSVertexerPPZ::EvalZ(TH1F *hist,Int_t sepa, Int_t ncoinc, TArrayF *zval) {
  fZFound=0;
  fZsig=0;
  Int_t N=0;
  Int_t NbinNotZero=0;
  Float_t totst = 0.;
  Float_t totst2 = 0.;
  Float_t curz = 0.;
  for(Int_t i=1;i<=sepa;i++){
    Float_t cont=hist->GetBinContent(i);
    curz = hist->GetBinLowEdge(i)+0.5*hist->GetBinWidth(i);
    totst+=cont;
    totst2+=cont*cont;
    N++;
    if(cont!=0)NbinNotZero++;
  }
  if(NbinNotZero==0){fZFound=-100; fZsig=-100; return;}
  if(NbinNotZero==1){fZFound = curz; fZsig=-100; return;}
  totst2=TMath::Sqrt(totst2/(N-1)-totst*totst/N/(N-1));
  totst /= N;
  Float_t val1=hist->GetBinLowEdge(sepa); 
  Float_t val2=hist->GetBinLowEdge(1);
  for(Int_t i=1;i<=sepa;i++){
    Float_t cont=hist->GetBinContent(i);
    if(cont>(totst+totst2*2.)){
      curz=hist->GetBinLowEdge(i);
      if(curz<val1)val1=curz;
      if(curz>val2)val2=curz;
    }
  }
  val2+=hist->GetBinWidth(1);
  if(fDebug>0)cout<<"Values for Z finding: "<<val1<<" "<<val2<<endl;
  fZFound=0;
  fZsig=0;
  N=0;
  for(Int_t i=0; i<ncoinc; i++){
    Float_t z=zval->At(i);
    if(z<val1)continue;
    if(z>val2)continue;
    fZFound+=z;
    fZsig+=z*z;
    N++;
  }
  if(N<=1){fZFound=-110; fZsig=-110; return;}
  fZsig=TMath::Sqrt((fZsig/(N-1)-fZFound*fZFound/N/(N-1))/N);
  fZFound=fZFound/N;
  fCurrentVertex = new AliITSVertex(fZFound,fZsig,N);
}

//______________________________________________________________________
AliITSVertex* AliITSVertexerPPZ::FindVertexForCurrentEvent(Int_t evnumber){
  // Defines the AliITSVertex for the current event
  fCurrentVertex = 0;
  fZFound = -999;
  fZsig = -999;
  if(!gAlice)gAlice = (AliRun*)fInFile->Get("gAlice");
  if(!gAlice){
    Error("FindVertexForCurrentEvent","The AliRun object is not available - nothing done");
    return fCurrentVertex;
  }
  if(!fITS)  {
    fITS = (AliITS*)gAlice->GetModule("ITS");
    if(!fITS) {
      Error("FindVertexForCurrentEvent","AliITS object was not found");
      return fCurrentVertex;
    }
  }
  AliITSgeom *geom = fITS->GetITSgeom();
  if(!geom) {
    Error("FindVertexForCurrentEvent","ITS geometry is not defined");
    return fCurrentVertex;
  }
  TTree *TR=0;
  TClonesArray *ITSrec  = 0;
  Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
  Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
  Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
  Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;

  TR = gAlice->TreeR();
  if(!TR){
    Error("FindVertexForCurrentEvent","TreeR not found");
    return fCurrentVertex;
  }
  ITSrec  = fITS->RecPoints();  // uses slow points or fast points if slow are
  // missing
  TBranch *branch = TR->GetBranch("ITSRecPoints");
  if(!branch){ 
    branch = TR->GetBranch("ITSRecPointsF");
    //   Warning("FindVertexForCurrentEvent","Using Fast points");
  }
  if(!branch){
   Error("FindVertexForCurrentEvent","branch for ITS rec points not found");
   return fCurrentVertex;
  }
  Float_t zave=0;
  Float_t rmszav=0;
  Float_t zave2=0;
  Int_t firipixe=0;
  for(Int_t module= fFirstL1; module<=fLastL1;module++){
    branch->GetEvent(module);
    Int_t nrecp1 = ITSrec->GetEntries();
    for(Int_t i=0; i<nrecp1;i++){
      AliITSRecPoint *current = (AliITSRecPoint*)ITSrec->At(i);
      lc[0]=current->GetX();
      lc[2]=current->GetZ();
      geom->LtoG(module,lc,gc);
      zave+=gc[2];
      zave2+=gc[2]*gc[2];
      firipixe++;
    }
    fITS->ResetRecPoints();
  }
  if(firipixe>1){
    rmszav=TMath::Sqrt(zave2/(firipixe-1)-zave*zave/firipixe/(firipixe-1));
    zave=zave/firipixe;
  }
  else {
    fZFound = -200;
    fZsig = -200;
    Warning("FindVertexForCurrentEvent","No rec points on first layer for this event");
    return fCurrentVertex;
  }
  Float_t zlim1=zave-rmszav;
  Float_t zlim2=zave+rmszav;
  Int_t sepa=(Int_t)((zlim2-zlim1)*10.+1.);
  zlim2=zlim1 + sepa/10.;
  TH1F *zvdis = new TH1F("z_ev","zv distr",sepa,zlim1,zlim2);
  if(fDebug>0)cout<<"Z limits: "<<zlim1<<" "<<zlim2<<"; Bins= "<<sepa<<endl;
  Int_t sizarr=100;
  TArrayF *zval = new TArrayF(sizarr);
  Int_t ncoinc=0;
  for(Int_t module= fFirstL1; module<=fLastL1;module++){
    if(fDebug>0)cout<<"processing module   "<<module<<"                  \r";
    branch->GetEvent(module);
    Int_t nrecp1 = ITSrec->GetEntries();
    TObjArray *poiL1 = new TObjArray(nrecp1);
    for(Int_t i=0; i<nrecp1;i++)poiL1->AddAt(ITSrec->At(i),i);
    fITS->ResetRecPoints();
    for(Int_t i=0; i<nrecp1;i++){
      AliITSRecPoint *current = (AliITSRecPoint*)poiL1->At(i);
      lc[0]=current->GetX();
      lc[2]=current->GetZ();
      geom->LtoG(module,lc,gc);
      gc[0]-=fX0;
      gc[1]-=fY0;
      Float_t r1=TMath::Sqrt(gc[0]*gc[0]+gc[1]*gc[1]);
      Float_t phi1 = TMath::ATan2(gc[1],gc[0]);
      if(phi1<0)phi1=2*TMath::Pi()+phi1;
      if(fDebug>1)cout<<"module "<<module<<" "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" "<<phi1<<"     \n";
      for(Int_t modul2=fFirstL2; modul2<=fLastL2; modul2++){
	branch->GetEvent(modul2);
	Int_t nrecp2 = ITSrec->GetEntries();
	for(Int_t j=0; j<nrecp2;j++){
	  AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(j);
	  lc2[0]=recp->GetX();
	  lc2[2]=recp->GetZ();
	  geom->LtoG(modul2,lc2,gc2);
	  gc2[0]-=fX0;
	  gc2[1]-=fY0;
	  Float_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
	  Float_t zr0=(r2*gc[2]-r1*gc2[2])/(r2-r1);
	  Float_t phi2 = TMath::ATan2(gc2[1],gc2[0]);
	  if(phi2<0)phi2=2.*TMath::Pi()+phi2;
	  Float_t diff = TMath::Abs(phi2-phi1);
	  if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
	  if(zr0>zlim1 && zr0<zlim2){
	    if(diff<fDiffPhiMax ){
	      zvdis->Fill(zr0);
	      zval->AddAt(zr0,ncoinc);
	      ncoinc++;
	      if(ncoinc==(sizarr-1)){
		sizarr+=100;
		zval->Set(sizarr);
	      }
	    }
	  }
	}
	fITS->ResetRecPoints();
      }
    }
    delete poiL1;
  }         // loop on modules
  if(fDebug>0){
    cout<<endl<<"Number of coincidences = "<<ncoinc<<endl;
  }
  EvalZ(zvdis,sepa,ncoinc,zval);
  delete zvdis;
  delete zval;
  if(fCurrentVertex){
    char name[30];
    sprintf(name,"Vertex_%d",evnumber);
    fCurrentVertex->SetName(name);
    fCurrentVertex->SetTitle("vertexer: PPZ");
  }
  return fCurrentVertex;
}

//______________________________________________________________________
void AliITSVertexerPPZ::FindVertices(){
  // computes the vertices of the events in the range FirstEvent - LastEvent
  if(!fOutFile){
    Error("FindVertices","The output file is not available - nothing done");
    return;
  }
  if(!gAlice)gAlice = (AliRun*)fInFile->Get("gAlice");
  if(!gAlice){
    Error("FindVertices","The AliRun object is not available - nothing done");
    return;
  }
  TDirectory *curdir = gDirectory;
  fInFile->cd();
  for(Int_t i=fFirstEvent;i<fLastEvent;i++){
    gAlice->GetEvent(i);
    FindVertexForCurrentEvent(i);
    if(fCurrentVertex){
      WriteCurrentVertex();
    }
    else {
      if(fDebug>0){
	cout<<"Vertex not found for event "<<i<<endl;
	cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
      }
    }
  }
  curdir->cd();
}

//________________________________________________________
void AliITSVertexerPPZ::PrintStatus() const {
  // Print current status
  cout <<"=======================================================\n";
  cout <<" First layer first and last modules: "<<fFirstL1<<", ";
  cout <<fLastL1<<endl;
  cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
  cout <<fLastL2<<endl;
  cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
  cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
  cout <<" Debug flag: "<<fDebug<<endl;
  if(fInFile)cout <<" The input file is open\n";
  if(!fInFile)cout <<"The input file is not open\n";
  if(fOutFile)cout <<"The output file is open\n";
  if(fOutFile)cout <<"The output file not is open\n";
  cout<<"First event to be processed "<<fFirstEvent;
  cout<<"\n Last event to be processed "<<fLastEvent<<endl;
}
Commit	Line	Data
c5f0f3c1	1	#include <AliITSVertexerPPZ.h>
	2	#include <Riostream.h>
	3	#include <TArrayF.h>
	4	#include <TDirectory.h>
	5	#include <TH1F.h>
	6	#include <TMath.h>
	7	#include <TObjArray.h>
	8	#include <TTree.h>
	9	#include <TClonesArray.h>
	10	#include "AliRun.h"
	11	#include "AliITS.h"
	12	#include "AliITSgeom.h"
	13	#include "AliITSRecPoint.h"
	14	/////////////////////////////////////////////////////////////////////////
	15	// //
	16	// This class is intended to compute the Z coordinate of the //
	17	// primary vertex for p-p interactions, or in general when the //
	18	// number of secondaries is too slow to use the class //
	19	// AliITSVertexerIons, which in turn is optimized for A-A collsions //
	20	// Origin: masera@to.infn.it 9/12/2002 //
	21	// //
	22	/////////////////////////////////////////////////////////////////////////
	23	ClassImp(AliITSVertexerPPZ)
	24
	25
	26
	27	//______________________________________________________________________
	28	AliITSVertexerPPZ::AliITSVertexerPPZ():AliITSVertexer() {
	29	// Default Constructor
	30
	31	SetDiffPhiMax(0);
	32	fX0 = 0.;
	33	fY0 = 0.;
	34	SetFirstLayerModules(0);
	35	SetSecondLayerModules(0);
	36	fZFound = 0;
	37	fZsig = 0.;
	38	fITS = 0;
	39	}
	40
	41	AliITSVertexerPPZ::AliITSVertexerPPZ(TFile infile, TFile outfile, Float_t x0, Float_t y0):AliITSVertexer(infile,outfile) {
	42	// Standard constructor
	43	if(!fInFile)Fatal("AliITSVertexerPPZ","No inputfile has been defined!");
	44	SetDiffPhiMax();
	45	fX0 = x0;
	46	fY0 = y0;
	47	SetFirstLayerModules();
	48	SetSecondLayerModules();
	49	fZFound = 0;
	50	fZsig = 0.;
	51	fITS = 0;
	52	}
	53
	54	//______________________________________________________________________
	55	AliITSVertexerPPZ::~AliITSVertexerPPZ() {
	56	// Default Destructor
	57	fITS = 0;
	58	}
	59
	60	//________________________________________________________
	61	void AliITSVertexerPPZ::EvalZ(TH1F hist,Int_t sepa, Int_t ncoinc, TArrayF zval) {
	62	fZFound=0;
	63	fZsig=0;
	64	Int_t N=0;
65	Int_t NbinNotZero=0;
66	Float_t totst = 0.;
67	Float_t totst2 = 0.;
68	Float_t curz = 0.;
69	for(Int_t i=1;i<=sepa;i++){
70	Float_t cont=hist->GetBinContent(i);
71	curz = hist->GetBinLowEdge(i)+0.5*hist->GetBinWidth(i);
72	totst+=cont;
73	totst2+=cont*cont;
74	N++;
75	if(cont!=0)NbinNotZero++;
76	}
77	if(NbinNotZero==0){fZFound=-100; fZsig=-100; return;}
78	if(NbinNotZero==1){fZFound = curz; fZsig=-100; return;}
79	totst2=TMath::Sqrt(totst2/(N-1)-totst*totst/N/(N-1));
80	totst /= N;
81	Float_t val1=hist->GetBinLowEdge(sepa);
82	Float_t val2=hist->GetBinLowEdge(1);
83	for(Int_t i=1;i<=sepa;i++){
84	Float_t cont=hist->GetBinContent(i);
85	if(cont>(totst+totst2*2.)){
86	curz=hist->GetBinLowEdge(i);
87	if(curz<val1)val1=curz;
88	if(curz>val2)val2=curz;
89	}
90	}
91	val2+=hist->GetBinWidth(1);
92	if(fDebug>0)cout<<"Values for Z finding: "<<val1<<" "<<val2<<endl;
93	fZFound=0;
94	fZsig=0;
95	N=0;
96	for(Int_t i=0; i<ncoinc; i++){
97	Float_t z=zval->At(i);
98	if(z<val1)continue;
99	if(z>val2)continue;
100	fZFound+=z;
101	fZsig+=z*z;
102	N++;
103	}
104	if(N<=1){fZFound=-110; fZsig=-110; return;}
105	fZsig=TMath::Sqrt((fZsig/(N-1)-fZFound*fZFound/N/(N-1))/N);
106	fZFound=fZFound/N;
107	fCurrentVertex = new AliITSVertex(fZFound,fZsig,N);
108	}
109
110	//______________________________________________________________________
111	AliITSVertex* AliITSVertexerPPZ::FindVertexForCurrentEvent(Int_t evnumber){
112	// Defines the AliITSVertex for the current event
113	fCurrentVertex = 0;
114	fZFound = -999;
115	fZsig = -999;
116	if(!gAlice)gAlice = (AliRun*)fInFile->Get("gAlice");
117	if(!gAlice){
118	Error("FindVertexForCurrentEvent","The AliRun object is not available - nothing done");
119	return fCurrentVertex;
120	}
121	if(!fITS) {
122	fITS = (AliITS*)gAlice->GetModule("ITS");
123	if(!fITS) {
124	Error("FindVertexForCurrentEvent","AliITS object was not found");
125	return fCurrentVertex;
126	}
127	}
128	AliITSgeom *geom = fITS->GetITSgeom();
129	if(!geom) {
130	Error("FindVertexForCurrentEvent","ITS geometry is not defined");
131	return fCurrentVertex;
132	}
133	TTree *TR=0;
134	TClonesArray *ITSrec = 0;
135	Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
136	Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
137	Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
138	Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
139
140	TR = gAlice->TreeR();
141	if(!TR){
142	Error("FindVertexForCurrentEvent","TreeR not found");
143	return fCurrentVertex;
144	}
145	ITSrec = fITS->RecPoints(); // uses slow points or fast points if slow are
146	// missing
147	TBranch *branch = TR->GetBranch("ITSRecPoints");
148	if(!branch){
149	branch = TR->GetBranch("ITSRecPointsF");
150	// Warning("FindVertexForCurrentEvent","Using Fast points");
151	}
152	if(!branch){
153	Error("FindVertexForCurrentEvent","branch for ITS rec points not found");
154	return fCurrentVertex;
155	}
156	Float_t zave=0;
157	Float_t rmszav=0;
158	Float_t zave2=0;
159	Int_t firipixe=0;
160	for(Int_t module= fFirstL1; module<=fLastL1;module++){
161	branch->GetEvent(module);
162	Int_t nrecp1 = ITSrec->GetEntries();
163	for(Int_t i=0; i<nrecp1;i++){
164	AliITSRecPoint current = (AliITSRecPoint)ITSrec->At(i);
165	lc[0]=current->GetX();
166	lc[2]=current->GetZ();
167	geom->LtoG(module,lc,gc);
168	zave+=gc[2];
169	zave2+=gc[2]*gc[2];
170	firipixe++;
171	}
172	fITS->ResetRecPoints();
173	}
174	if(firipixe>1){
175	rmszav=TMath::Sqrt(zave2/(firipixe-1)-zave*zave/firipixe/(firipixe-1));
176	zave=zave/firipixe;
177	}
178	else {
179	fZFound = -200;
180	fZsig = -200;
181	Warning("FindVertexForCurrentEvent","No rec points on first layer for this event");
182	return fCurrentVertex;
183	}
184	Float_t zlim1=zave-rmszav;
185	Float_t zlim2=zave+rmszav;
186	Int_t sepa=(Int_t)((zlim2-zlim1)*10.+1.);
187	zlim2=zlim1 + sepa/10.;
188	TH1F *zvdis = new TH1F("z_ev","zv distr",sepa,zlim1,zlim2);
189	if(fDebug>0)cout<<"Z limits: "<<zlim1<<" "<<zlim2<<"; Bins= "<<sepa<<endl;
190	Int_t sizarr=100;
191	TArrayF *zval = new TArrayF(sizarr);
192	Int_t ncoinc=0;
193	for(Int_t module= fFirstL1; module<=fLastL1;module++){
194	if(fDebug>0)cout<<"processing module "<<module<<" \r";
195	branch->GetEvent(module);
196	Int_t nrecp1 = ITSrec->GetEntries();
197	TObjArray *poiL1 = new TObjArray(nrecp1);
198	for(Int_t i=0; i<nrecp1;i++)poiL1->AddAt(ITSrec->At(i),i);
199	fITS->ResetRecPoints();
200	for(Int_t i=0; i<nrecp1;i++){
201	AliITSRecPoint current = (AliITSRecPoint)poiL1->At(i);
202	lc[0]=current->GetX();
203	lc[2]=current->GetZ();
204	geom->LtoG(module,lc,gc);
205	gc[0]-=fX0;
206	gc[1]-=fY0;
207	Float_t r1=TMath::Sqrt(gc[0]gc[0]+gc[1]gc[1]);
208	Float_t phi1 = TMath::ATan2(gc[1],gc[0]);
209	if(phi1<0)phi1=2*TMath::Pi()+phi1;
210	if(fDebug>1)cout<<"module "<<module<<" "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" "<<phi1<<" \n";
211	for(Int_t modul2=fFirstL2; modul2<=fLastL2; modul2++){
212	branch->GetEvent(modul2);
213	Int_t nrecp2 = ITSrec->GetEntries();
214	for(Int_t j=0; j<nrecp2;j++){
215	AliITSRecPoint recp = (AliITSRecPoint)ITSrec->At(j);
216	lc2[0]=recp->GetX();
217	lc2[2]=recp->GetZ();
218	geom->LtoG(modul2,lc2,gc2);
219	gc2[0]-=fX0;
220	gc2[1]-=fY0;
221	Float_t r2=TMath::Sqrt(gc2[0]gc2[0]+gc2[1]gc2[1]);
222	Float_t zr0=(r2gc[2]-r1gc2[2])/(r2-r1);
223	Float_t phi2 = TMath::ATan2(gc2[1],gc2[0]);
224	if(phi2<0)phi2=2.*TMath::Pi()+phi2;
225	Float_t diff = TMath::Abs(phi2-phi1);
226	if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
227	if(zr0>zlim1 && zr0<zlim2){
228	if(diff<fDiffPhiMax ){
229	zvdis->Fill(zr0);
230	zval->AddAt(zr0,ncoinc);
231	ncoinc++;
232	if(ncoinc==(sizarr-1)){
233	sizarr+=100;
234	zval->Set(sizarr);
235	}
236	}
237	}
238	}
239	fITS->ResetRecPoints();
240	}
241	}
242	delete poiL1;
243	} // loop on modules
244	if(fDebug>0){
245	cout<<endl<<"Number of coincidences = "<<ncoinc<<endl;
246	}
247	EvalZ(zvdis,sepa,ncoinc,zval);
248	delete zvdis;
249	delete zval;
250	if(fCurrentVertex){
251	char name[30];
252	sprintf(name,"Vertex_%d",evnumber);
253	fCurrentVertex->SetName(name);
254	fCurrentVertex->SetTitle("vertexer: PPZ");
255	}
256	return fCurrentVertex;
257	}
258
259	//______________________________________________________________________
260	void AliITSVertexerPPZ::FindVertices(){
261	// computes the vertices of the events in the range FirstEvent - LastEvent
262	if(!fOutFile){
263	Error("FindVertices","The output file is not available - nothing done");
264	return;
265	}
266	if(!gAlice)gAlice = (AliRun*)fInFile->Get("gAlice");
267	if(!gAlice){
268	Error("FindVertices","The AliRun object is not available - nothing done");
269	return;
270	}
271	TDirectory *curdir = gDirectory;
272	fInFile->cd();
273	for(Int_t i=fFirstEvent;i<fLastEvent;i++){
274	gAlice->GetEvent(i);
275	FindVertexForCurrentEvent(i);
276	if(fCurrentVertex){
277	WriteCurrentVertex();
278	}
279	else {
280	if(fDebug>0){
281	cout<<"Vertex not found for event "<<i<<endl;
282	cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
283	}
284	}
285	}
286	curdir->cd();
287	}
288
289	//________________________________________________________
290	void AliITSVertexerPPZ::PrintStatus() const {
291	// Print current status
292	cout <<"=======================================================\n";
293	cout <<" First layer first and last modules: "<<fFirstL1<<", ";
294	cout <<fLastL1<<endl;
295	cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
296	cout <<fLastL2<<endl;
297	cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
298	cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
299	cout <<" Debug flag: "<<fDebug<<endl;
300	if(fInFile)cout <<" The input file is open\n";
301	if(!fInFile)cout <<"The input file is not open\n";
302	if(fOutFile)cout <<"The output file is open\n";
303	if(fOutFile)cout <<"The output file not is open\n";
304	cout<<"First event to be processed "<<fFirstEvent;
305	cout<<"\n Last event to be processed "<<fLastEvent<<endl;
306	}