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

/**************************************************************************
 * Copyright(c) 1998-2003, 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 <AliITSVertexerPPZ.h>
#include <TArrayF.h>
#include <TH1F.h>
#include <TObjArray.h>
#include <TTree.h>
#include <TClonesArray.h>
#include "AliITSDetTypeRec.h"
#include "AliITSgeom.h"
#include "AliITSLoader.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;
  SetWindow(0);
}

AliITSVertexerPPZ::AliITSVertexerPPZ(TString fn, Float_t x0, Float_t y0):AliITSVertexer(fn) {
  // Standard constructor
  SetDiffPhiMax();
  fX0 = x0;
  fY0 = y0;
  SetFirstLayerModules();
  SetSecondLayerModules();
  fZFound = 0;
  fZsig = 0.;
  //fITS = 0;
  SetWindow();
}

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

//________________________________________________________
void  AliITSVertexerPPZ::EvalZ(TH1F *hist,Int_t sepa, Int_t ncoinc, TArrayF *zval) {

  //Evaluation of Z
  Float_t deltaVal = hist->GetBinWidth(1)*fWindow; // max window in Z for searching
  fZFound=0;
  fZsig=0;
  Int_t nN=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);
    if(cont!=0)curz = hist->GetBinLowEdge(i)+0.5*hist->GetBinWidth(i);
    totst+=cont;
    totst2+=cont*cont;
    nN++;
    if(cont!=0)nbinNotZero++;
  }
  if(nbinNotZero==0){fZFound=-100; fZsig=-100; return;}
  if(nbinNotZero==1){
    fZFound = curz;
    fZsig=0;
    fCurrentVertex = new AliESDVertex(fZFound,fZsig,nbinNotZero);
    return;
  }
  Float_t errsq = totst2/(nN-1)-totst*totst/nN/(nN-1);
  if(errsq>=0){
  totst2=TMath::Sqrt(errsq);
  }
  else {
    Error("EvalZ","Negative variance: %d - %12.7f - %12.7f",nN,totst2,totst);
    fZFound=-100; 
    fZsig=-100;
    return;
  }
  totst /= nN;
  Float_t cut = totst+totst2*2.;
  if(fDebug>1)cout<<"totst, totst2, cut: "<<totst<<", "<<totst2<<", "<<cut<<endl;
  Float_t val1=hist->GetBinLowEdge(sepa); 
  Float_t val2=hist->GetBinLowEdge(1);
  Float_t valm = 0.;
  Float_t zmax = 0.;
  for(Int_t i=1;i<=sepa;i++){
    Float_t cont=hist->GetBinContent(i);
    if(cont>valm){
      valm = cont;
      zmax = hist->GetBinLowEdge(i)+0.5*hist->GetBinWidth(1);
    }
    if(cont>cut){
      curz=hist->GetBinLowEdge(i);
      if(curz<val1)val1=curz;
      if(curz>val2)val2=curz;
    }
  }
  val2+=hist->GetBinWidth(1);
  if((val2-val1)>deltaVal){
    val1 = zmax-deltaVal/2.;
    val2 = zmax+deltaVal/2.;
    if(fDebug>0)cout<<"val1 and val2 recomputed\n";
  }
  if(fDebug>0)cout<<"Values for Z finding: "<<val1<<" "<<val2<<" "<<val2-val1<<endl;
  fZFound=0;
  fZsig=0;
  nN=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;
    /*   weights defined by the curvature
    Float_t wei = 1./curv->At(i);
    fZFound+=z*wei;
    fZsig+=wei;
    */
    nN++;
  }
  if(nN<1){fZFound=-110; fZsig=-110; return;}
  if(nN==1){
    fZsig = 0;
    fCurrentVertex = new AliESDVertex(fZFound,fZsig,nN);
    return;
  }
  errsq = (fZsig/(nN-1)-fZFound*fZFound/nN/(nN-1))/nN;
  if(errsq>=0.){
    fZsig=TMath::Sqrt(errsq);
  }
  else {
    Error("evalZ","Negative variance: %d - %12.7f %12.7f",nN,fZsig,fZFound);
    fZsig=0;
  }
  fZFound=fZFound/nN;
  /* weights defined by the curvature
  fZsig=1./fZsig;
  fZFound*=fZsig;
  fZsig = TMath::Sqrt(fZsig);
  */
  fCurrentVertex = new AliESDVertex(fZFound,fZsig,nN);
}

//______________________________________________________________________
AliESDVertex* AliITSVertexerPPZ::FindVertexForCurrentEvent(Int_t evnumber){
  // Defines the AliESDVertex for the current event
  fCurrentVertex = 0;
  fZFound = -999;
  fZsig = -999;
  AliRunLoader *rl =AliRunLoader::GetRunLoader();
  AliITSLoader* iTSloader = (AliITSLoader*)rl->GetLoader("ITSLoader");
  TDirectory * olddir = gDirectory;
  rl->CdGAFile();
  AliITSgeom* geom = (AliITSgeom*)gDirectory->Get("AliITSgeom");
  olddir->cd(); 

  AliITSDetTypeRec detTypeRec;

  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 = iTSloader->TreeR();
  if(!tR){
    Error("FindVertexForCurrentEvent","TreeR not found");
    return fCurrentVertex;
  }
  detTypeRec.SetTreeAddressR(tR);
  itsRec = detTypeRec.RecPoints();
  // missing
  // TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
  TBranch *branch;
  //if(fUseV2Clusters){
  //  branch = tR->GetBranch("ITSRecPoints");
  //  branch->SetAddress(&clusters);
  //}
  //else {
  branch = tR->GetBranch("ITSRecPoints");
  if(!branch){ 
    branch = tR->GetBranch("ITSRecPointsF");
  }
  //}
  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);
    //if(fUseV2Clusters){
    //  Clusters2RecPoints(clusters,module,itsRec);
    //}
    Int_t nrecp1 = itsRec->GetEntries();
    for(Int_t i=0; i<nrecp1;i++){
      AliITSRecPoint *current = (AliITSRecPoint*)itsRec->At(i);
      lc[0]=current->GetDetLocalX();
      lc[2]=current->GetDetLocalZ();
      geom->LtoG(module,lc,gc);
      zave+=gc[2];
      zave2+=gc[2]*gc[2];
      firipixe++;
    }
    detTypeRec.ResetRecPoints();
  }
  if(firipixe>1){
    rmszav=TMath::Sqrt(zave2/(firipixe-1)-zave*zave/firipixe/(firipixe-1));
    zave=zave/firipixe;
    if(fDebug>1)cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++\n Number of firing pixels: "<<firipixe<<endl;
  }
  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;
    cout<<"Bin width: "<<zvdis->GetBinWidth(1)*10000.<<" micron\n";
  }
  Int_t sizarr=100;
  TArrayF *zval = new TArrayF(sizarr);
  //  TArrayF *curv = 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);
    //if(fUseV2Clusters){
    //  Clusters2RecPoints(clusters,module,itsRec);
    //}
    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();
    detTypeRec.ResetRecPoints();
    for(Int_t i=0; i<nrecp1;i++){
      AliITSRecPoint *current = (AliITSRecPoint*)poiL1->At(i);
      lc[0]=current->GetDetLocalX();
      lc[2]=current->GetDetLocalZ();
      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);
	//if(fUseV2Clusters){
	//  Clusters2RecPoints(clusters,modul2,itsRec);
	//}
	Int_t nrecp2 = itsRec->GetEntries();
	for(Int_t j=0; j<nrecp2;j++){
	  AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
	  lc2[0]=recp->GetDetLocalX();
	  lc2[2]=recp->GetDetLocalZ();
	  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);
	      /* uncomment these lines to use curvature as a weight
	      Float_t cu = Curv(0.,0.,gc[0],gc[1],gc2[0],gc[1]);
	      curv->AddAt(cu,ncoinc);
	      */
	      ncoinc++;
	      if(ncoinc==(sizarr-1)){
		sizarr+=100;
		zval->Set(sizarr);
		//uncomment next line to use curvature as weight
		//		curv->Set(sizarr);
	      }
	    }
	  }
	}
	detTypeRec.ResetRecPoints();
      }
    }
    delete poiL1;
  }         // loop on modules
  if(fDebug>0){
    cout<<endl<<"Number of coincidences = "<<ncoinc<<endl;
  }
  //  EvalZ(zvdis,sepa,ncoinc,zval,curv);
  EvalZ(zvdis,sepa,ncoinc,zval);
  delete zvdis;
  delete zval;
  //  delete curv;
  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
  AliRunLoader *rl = AliRunLoader::GetRunLoader();
  AliITSLoader* iTSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
  iTSloader->ReloadRecPoints();
  for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
    cout<<"Processing event "<<i<<endl;
    rl->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;
      }
    }
  }
}

//________________________________________________________
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 <<" Window for Z search: "<<fWindow<<endl;
  cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
  cout <<" Debug flag: "<<fDebug<<endl;
  cout<<"First event to be processed "<<fFirstEvent;
  cout<<"\n Last event to be processed "<<fLastEvent<<endl;
}

//________________________________________________________
Float_t AliITSVertexerPPZ::Curv(Double_t x1,Double_t y1,
                   Double_t x2,Double_t y2,
                   Double_t x3,Double_t y3) 
{
  //-----------------------------------------------------------------
  // Initial approximation of the track curvature (Y. Belikov) squared
  //-----------------------------------------------------------------
  Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
  Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
                  (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
  Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
                  (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));

  Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);

  Float_t val = (Float_t)(( -xr*yr/TMath::Sqrt(xr*xr+yr*yr))
			  *( -xr*yr/TMath::Sqrt(xr*xr+yr*yr)));
  return val; 
}
Commit	Line	Data
c95d6417	1	/**************************************************************************
	2	* Copyright(c) 1998-2003, 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	**************************************************************************/
c5f0f3c1	15	#include <AliITSVertexerPPZ.h>
c5f0f3c1	16	#include <TArrayF.h>
c5f0f3c1	17	#include <TH1F.h>
c5f0f3c1	18	#include <TObjArray.h>
	19	#include <TTree.h>
	20	#include <TClonesArray.h>
7d62fb64	21	#include "AliITSDetTypeRec.h"
c5f0f3c1	22	#include "AliITSgeom.h"
88cb7938	23	#include "AliITSLoader.h"
c5f0f3c1	24	#include "AliITSRecPoint.h"
	25	/////////////////////////////////////////////////////////////////////////
	26	// //
	27	// This class is intended to compute the Z coordinate of the //
	28	// primary vertex for p-p interactions, or in general when the //
	29	// number of secondaries is too slow to use the class //
	30	// AliITSVertexerIons, which in turn is optimized for A-A collsions //
	31	// Origin: masera@to.infn.it 9/12/2002 //
	32	// //
	33	/////////////////////////////////////////////////////////////////////////
	34	ClassImp(AliITSVertexerPPZ)
	35
	36
	37
	38	//______________________________________________________________________
	39	AliITSVertexerPPZ::AliITSVertexerPPZ():AliITSVertexer() {
	40	// Default Constructor
	41
	42	SetDiffPhiMax(0);
	43	fX0 = 0.;
	44	fY0 = 0.;
	45	SetFirstLayerModules(0);
	46	SetSecondLayerModules(0);
	47	fZFound = 0;
	48	fZsig = 0.;
7d62fb64	49	//fITS = 0;
c95d6417	50	SetWindow(0);
c5f0f3c1	51	}
c5f0f3c1	52
88cb7938	53	AliITSVertexerPPZ::AliITSVertexerPPZ(TString fn, Float_t x0, Float_t y0):AliITSVertexer(fn) {
c5f0f3c1	54	// Standard constructor
c5f0f3c1	55	SetDiffPhiMax();
	56	fX0 = x0;
	57	fY0 = y0;
	58	SetFirstLayerModules();
	59	SetSecondLayerModules();
	60	fZFound = 0;
	61	fZsig = 0.;
7d62fb64	62	//fITS = 0;
c95d6417	63	SetWindow();
c5f0f3c1	64	}
	65
	66	//______________________________________________________________________
	67	AliITSVertexerPPZ::~AliITSVertexerPPZ() {
	68	// Default Destructor
7d62fb64	69	//fITS = 0;
c5f0f3c1	70	}
	71
	72	//________________________________________________________
	73	void AliITSVertexerPPZ::EvalZ(TH1F hist,Int_t sepa, Int_t ncoinc, TArrayF zval) {
7d62fb64	74
	75	//Evaluation of Z
	76	Float_t deltaVal = hist->GetBinWidth(1)*fWindow; // max window in Z for searching
c5f0f3c1	77	fZFound=0;
c5f0f3c1	78	fZsig=0;
7d62fb64	79	Int_t nN=0;
7d62fb64	80	Int_t nbinNotZero=0;
c5f0f3c1	81	Float_t totst = 0.;
	82	Float_t totst2 = 0.;
	83	Float_t curz = 0.;
	84	for(Int_t i=1;i<=sepa;i++){
	85	Float_t cont=hist->GetBinContent(i);
c95d6417	86	if(cont!=0)curz = hist->GetBinLowEdge(i)+0.5*hist->GetBinWidth(i);
c5f0f3c1	87	totst+=cont;
c5f0f3c1	88	totst2+=cont*cont;
7d62fb64	89	nN++;
7d62fb64	90	if(cont!=0)nbinNotZero++;
c5f0f3c1	91	}
7d62fb64	92	if(nbinNotZero==0){fZFound=-100; fZsig=-100; return;}
7d62fb64	93	if(nbinNotZero==1){
c95d6417	94	fZFound = curz;
c95d6417	95	fZsig=0;
7d62fb64	96	fCurrentVertex = new AliESDVertex(fZFound,fZsig,nbinNotZero);
c95d6417	97	return;
c95d6417	98	}
7d62fb64	99	Float_t errsq = totst2/(nN-1)-totst*totst/nN/(nN-1);
c95d6417	100	if(errsq>=0){
	101	totst2=TMath::Sqrt(errsq);
	102	}
	103	else {
7d62fb64	104	Error("EvalZ","Negative variance: %d - %12.7f - %12.7f",nN,totst2,totst);
c95d6417	105	fZFound=-100;
	106	fZsig=-100;
	107	return;
	108	}
7d62fb64	109	totst /= nN;
c95d6417	110	Float_t cut = totst+totst2*2.;
c95d6417	111	if(fDebug>1)cout<<"totst, totst2, cut: "<<totst<<", "<<totst2<<", "<<cut<<endl;
c5f0f3c1	112	Float_t val1=hist->GetBinLowEdge(sepa);
c5f0f3c1	113	Float_t val2=hist->GetBinLowEdge(1);
c95d6417	114	Float_t valm = 0.;
c95d6417	115	Float_t zmax = 0.;
c5f0f3c1	116	for(Int_t i=1;i<=sepa;i++){
c5f0f3c1	117	Float_t cont=hist->GetBinContent(i);
c95d6417	118	if(cont>valm){
	119	valm = cont;
	120	zmax = hist->GetBinLowEdge(i)+0.5*hist->GetBinWidth(1);
	121	}
	122	if(cont>cut){
c5f0f3c1	123	curz=hist->GetBinLowEdge(i);
	124	if(curz<val1)val1=curz;
	125	if(curz>val2)val2=curz;
	126	}
	127	}
	128	val2+=hist->GetBinWidth(1);
7d62fb64	129	if((val2-val1)>deltaVal){
	130	val1 = zmax-deltaVal/2.;
	131	val2 = zmax+deltaVal/2.;
c95d6417	132	if(fDebug>0)cout<<"val1 and val2 recomputed\n";
	133	}
	134	if(fDebug>0)cout<<"Values for Z finding: "<<val1<<" "<<val2<<" "<<val2-val1<<endl;
c5f0f3c1	135	fZFound=0;
c5f0f3c1	136	fZsig=0;
7d62fb64	137	nN=0;
c5f0f3c1	138	for(Int_t i=0; i<ncoinc; i++){
	139	Float_t z=zval->At(i);
	140	if(z<val1)continue;
	141	if(z>val2)continue;
c95d6417	142
c5f0f3c1	143	fZFound+=z;
c5f0f3c1	144	fZsig+=z*z;
c95d6417	145	/* weights defined by the curvature
	146	Float_t wei = 1./curv->At(i);
	147	fZFound+=z*wei;
	148	fZsig+=wei;
	149	*/
7d62fb64	150	nN++;
c5f0f3c1	151	}
7d62fb64	152	if(nN<1){fZFound=-110; fZsig=-110; return;}
7d62fb64	153	if(nN==1){
c95d6417	154	fZsig = 0;
7d62fb64	155	fCurrentVertex = new AliESDVertex(fZFound,fZsig,nN);
c95d6417	156	return;
c95d6417	157	}
7d62fb64	158	errsq = (fZsig/(nN-1)-fZFound*fZFound/nN/(nN-1))/nN;
c95d6417	159	if(errsq>=0.){
	160	fZsig=TMath::Sqrt(errsq);
	161	}
	162	else {
7d62fb64	163	Error("evalZ","Negative variance: %d - %12.7f %12.7f",nN,fZsig,fZFound);
c95d6417	164	fZsig=0;
c95d6417	165	}
7d62fb64	166	fZFound=fZFound/nN;
c95d6417	167	/* weights defined by the curvature
	168	fZsig=1./fZsig;
	169	fZFound*=fZsig;
	170	fZsig = TMath::Sqrt(fZsig);
	171	*/
7d62fb64	172	fCurrentVertex = new AliESDVertex(fZFound,fZsig,nN);
c5f0f3c1	173	}
	174
	175	//______________________________________________________________________
d681bb2d	176	AliESDVertex* AliITSVertexerPPZ::FindVertexForCurrentEvent(Int_t evnumber){
d681bb2d	177	// Defines the AliESDVertex for the current event
c5f0f3c1	178	fCurrentVertex = 0;
	179	fZFound = -999;
	180	fZsig = -999;
88cb7938	181	AliRunLoader *rl =AliRunLoader::GetRunLoader();
7d62fb64	182	AliITSLoader* iTSloader = (AliITSLoader*)rl->GetLoader("ITSLoader");
	183	TDirectory * olddir = gDirectory;
	184	rl->CdGAFile();
	185	AliITSgeom* geom = (AliITSgeom*)gDirectory->Get("AliITSgeom");
	186	olddir->cd();
	187
	188	AliITSDetTypeRec detTypeRec;
	189
c5f0f3c1	190	if(!geom) {
	191	Error("FindVertexForCurrentEvent","ITS geometry is not defined");
	192	return fCurrentVertex;
	193	}
2257f27e	194	TTree *tR=0;
2257f27e	195	TClonesArray *itsRec = 0;
c5f0f3c1	196	Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
	197	Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
	198	Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
	199	Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
	200
7d62fb64	201	tR = iTSloader->TreeR();
2257f27e	202	if(!tR){
c5f0f3c1	203	Error("FindVertexForCurrentEvent","TreeR not found");
	204	return fCurrentVertex;
	205	}
7d62fb64	206	detTypeRec.SetTreeAddressR(tR);
7d62fb64	207	itsRec = detTypeRec.RecPoints();
c5f0f3c1	208	// missing
00a7cc50	209	// TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
2257f27e	210	TBranch *branch;
00a7cc50	211	//if(fUseV2Clusters){
	212	// branch = tR->GetBranch("ITSRecPoints");
	213	// branch->SetAddress(&clusters);
	214	//}
	215	//else {
	216	branch = tR->GetBranch("ITSRecPoints");
	217	if(!branch){
	218	branch = tR->GetBranch("ITSRecPointsF");
c5f0f3c1	219	}
00a7cc50	220	//}
c5f0f3c1	221	if(!branch){
	222	Error("FindVertexForCurrentEvent","branch for ITS rec points not found");
	223	return fCurrentVertex;
	224	}
	225	Float_t zave=0;
	226	Float_t rmszav=0;
	227	Float_t zave2=0;
	228	Int_t firipixe=0;
	229	for(Int_t module= fFirstL1; module<=fLastL1;module++){
	230	branch->GetEvent(module);
00a7cc50	231	//if(fUseV2Clusters){
	232	// Clusters2RecPoints(clusters,module,itsRec);
	233	//}
2257f27e	234	Int_t nrecp1 = itsRec->GetEntries();
c5f0f3c1	235	for(Int_t i=0; i<nrecp1;i++){
2257f27e	236	AliITSRecPoint current = (AliITSRecPoint)itsRec->At(i);
00a7cc50	237	lc[0]=current->GetDetLocalX();
00a7cc50	238	lc[2]=current->GetDetLocalZ();
c5f0f3c1	239	geom->LtoG(module,lc,gc);
	240	zave+=gc[2];
	241	zave2+=gc[2]*gc[2];
	242	firipixe++;
	243	}
7d62fb64	244	detTypeRec.ResetRecPoints();
c5f0f3c1	245	}
	246	if(firipixe>1){
	247	rmszav=TMath::Sqrt(zave2/(firipixe-1)-zave*zave/firipixe/(firipixe-1));
	248	zave=zave/firipixe;
c95d6417	249	if(fDebug>1)cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++\n Number of firing pixels: "<<firipixe<<endl;
c5f0f3c1	250	}
	251	else {
	252	fZFound = -200;
	253	fZsig = -200;
	254	Warning("FindVertexForCurrentEvent","No rec points on first layer for this event");
	255	return fCurrentVertex;
	256	}
	257	Float_t zlim1=zave-rmszav;
	258	Float_t zlim2=zave+rmszav;
	259	Int_t sepa=(Int_t)((zlim2-zlim1)*10.+1.);
	260	zlim2=zlim1 + sepa/10.;
	261	TH1F *zvdis = new TH1F("z_ev","zv distr",sepa,zlim1,zlim2);
c95d6417	262	if(fDebug>0){
	263	cout<<"Z limits: "<<zlim1<<" "<<zlim2<<"; Bins= "<<sepa<<endl;
	264	cout<<"Bin width: "<<zvdis->GetBinWidth(1)*10000.<<" micron\n";
	265	}
c5f0f3c1	266	Int_t sizarr=100;
c5f0f3c1	267	TArrayF *zval = new TArrayF(sizarr);
c95d6417	268	// TArrayF *curv = new TArrayF(sizarr);
c5f0f3c1	269	Int_t ncoinc=0;
	270	for(Int_t module= fFirstL1; module<=fLastL1;module++){
	271	if(fDebug>0)cout<<"processing module "<<module<<" \r";
	272	branch->GetEvent(module);
00a7cc50	273	//if(fUseV2Clusters){
	274	// Clusters2RecPoints(clusters,module,itsRec);
	275	//}
2257f27e	276	Int_t nrecp1 = itsRec->GetEntries();
c5f0f3c1	277	TObjArray *poiL1 = new TObjArray(nrecp1);
2257f27e	278	for(Int_t i=0; i<nrecp1;i++)poiL1->AddAt(itsRec->At(i),i);
7d62fb64	279	//fITS->ResetRecPoints();
7d62fb64	280	detTypeRec.ResetRecPoints();
c5f0f3c1	281	for(Int_t i=0; i<nrecp1;i++){
c5f0f3c1	282	AliITSRecPoint current = (AliITSRecPoint)poiL1->At(i);
00a7cc50	283	lc[0]=current->GetDetLocalX();
00a7cc50	284	lc[2]=current->GetDetLocalZ();
c5f0f3c1	285	geom->LtoG(module,lc,gc);
	286	gc[0]-=fX0;
	287	gc[1]-=fY0;
	288	Float_t r1=TMath::Sqrt(gc[0]gc[0]+gc[1]gc[1]);
	289	Float_t phi1 = TMath::ATan2(gc[1],gc[0]);
	290	if(phi1<0)phi1=2*TMath::Pi()+phi1;
	291	if(fDebug>1)cout<<"module "<<module<<" "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" "<<phi1<<" \n";
	292	for(Int_t modul2=fFirstL2; modul2<=fLastL2; modul2++){
	293	branch->GetEvent(modul2);
00a7cc50	294	//if(fUseV2Clusters){
	295	// Clusters2RecPoints(clusters,modul2,itsRec);
	296	//}
2257f27e	297	Int_t nrecp2 = itsRec->GetEntries();
c5f0f3c1	298	for(Int_t j=0; j<nrecp2;j++){
2257f27e	299	AliITSRecPoint recp = (AliITSRecPoint)itsRec->At(j);
00a7cc50	300	lc2[0]=recp->GetDetLocalX();
00a7cc50	301	lc2[2]=recp->GetDetLocalZ();
c5f0f3c1	302	geom->LtoG(modul2,lc2,gc2);
	303	gc2[0]-=fX0;
	304	gc2[1]-=fY0;
	305	Float_t r2=TMath::Sqrt(gc2[0]gc2[0]+gc2[1]gc2[1]);
	306	Float_t zr0=(r2gc[2]-r1gc2[2])/(r2-r1);
	307	Float_t phi2 = TMath::ATan2(gc2[1],gc2[0]);
	308	if(phi2<0)phi2=2.*TMath::Pi()+phi2;
	309	Float_t diff = TMath::Abs(phi2-phi1);
	310	if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
	311	if(zr0>zlim1 && zr0<zlim2){
	312	if(diff<fDiffPhiMax ){
	313	zvdis->Fill(zr0);
	314	zval->AddAt(zr0,ncoinc);
c95d6417	315	/* uncomment these lines to use curvature as a weight
	316	Float_t cu = Curv(0.,0.,gc[0],gc[1],gc2[0],gc[1]);
	317	curv->AddAt(cu,ncoinc);
	318	*/
c5f0f3c1	319	ncoinc++;
	320	if(ncoinc==(sizarr-1)){
	321	sizarr+=100;
	322	zval->Set(sizarr);
c95d6417	323	//uncomment next line to use curvature as weight
c95d6417	324	// curv->Set(sizarr);
c5f0f3c1	325	}
	326	}
	327	}
	328	}
7d62fb64	329	detTypeRec.ResetRecPoints();
c5f0f3c1	330	}
	331	}
	332	delete poiL1;
	333	} // loop on modules
	334	if(fDebug>0){
	335	cout<<endl<<"Number of coincidences = "<<ncoinc<<endl;
	336	}
c95d6417	337	// EvalZ(zvdis,sepa,ncoinc,zval,curv);
c5f0f3c1	338	EvalZ(zvdis,sepa,ncoinc,zval);
	339	delete zvdis;
	340	delete zval;
c95d6417	341	// delete curv;
c5f0f3c1	342	if(fCurrentVertex){
	343	char name[30];
	344	sprintf(name,"Vertex_%d",evnumber);
88cb7938	345	// fCurrentVertex->SetName(name);
c5f0f3c1	346	fCurrentVertex->SetTitle("vertexer: PPZ");
c5f0f3c1	347	}
7d62fb64	348
c5f0f3c1	349	return fCurrentVertex;
	350	}
	351
	352	//______________________________________________________________________
	353	void AliITSVertexerPPZ::FindVertices(){
	354	// computes the vertices of the events in the range FirstEvent - LastEvent
88cb7938	355	AliRunLoader *rl = AliRunLoader::GetRunLoader();
7d62fb64	356	AliITSLoader* iTSloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
7d62fb64	357	iTSloader->ReloadRecPoints();
ab6a6f40	358	for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
88cb7938	359	cout<<"Processing event "<<i<<endl;
88cb7938	360	rl->GetEvent(i);
c5f0f3c1	361	FindVertexForCurrentEvent(i);
	362	if(fCurrentVertex){
	363	WriteCurrentVertex();
	364	}
	365	else {
	366	if(fDebug>0){
	367	cout<<"Vertex not found for event "<<i<<endl;
	368	cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
	369	}
	370	}
	371	}
c5f0f3c1	372	}
	373
	374	//________________________________________________________
	375	void AliITSVertexerPPZ::PrintStatus() const {
	376	// Print current status
	377	cout <<"=======================================================\n";
	378	cout <<" First layer first and last modules: "<<fFirstL1<<", ";
	379	cout <<fLastL1<<endl;
	380	cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
	381	cout <<fLastL2<<endl;
	382	cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
c95d6417	383	cout <<" Window for Z search: "<<fWindow<<endl;
c5f0f3c1	384	cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
c5f0f3c1	385	cout <<" Debug flag: "<<fDebug<<endl;
c5f0f3c1	386	cout<<"First event to be processed "<<fFirstEvent;
	387	cout<<"\n Last event to be processed "<<fLastEvent<<endl;
	388	}
c95d6417	389
	390	//________________________________________________________
	391	Float_t AliITSVertexerPPZ::Curv(Double_t x1,Double_t y1,
	392	Double_t x2,Double_t y2,
	393	Double_t x3,Double_t y3)
	394	{
	395	//-----------------------------------------------------------------
	396	// Initial approximation of the track curvature (Y. Belikov) squared
	397	//-----------------------------------------------------------------
	398	Double_t d=(x2-x1)(y3-y2)-(x3-x2)(y2-y1);
	399	Double_t a=0.5((y3-y2)(y2y2-y1y1+x2x2-x1x1)-
	400	(y2-y1)(y3y3-y2y2+x3x3-x2*x2));
	401	Double_t b=0.5((x2-x1)(y3y3-y2y2+x3x3-x2x2)-
	402	(x3-x2)(y2y2-y1y1+x2x2-x1*x1));
	403
	404	Double_t xr=TMath::Abs(d/(dx1-a)), yr=d/(dy1-b);
	405
	406	Float_t val = (Float_t)(( -xryr/TMath::Sqrt(xrxr+yr*yr))
	407	( -xryr/TMath::Sqrt(xrxr+yryr)));
	408	return val;
	409	}