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

/**************************************************************************
 * Copyright(c) 2006-2008, 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.                  *
 **************************************************************************/

//-----------------------------------------------------------------
//           AliITSVertexer3DTapan class
//   This is a class for the 3d vertex finding
//    Origin: Tapan Nayak, VECC-CERN, Tapan.Nayak@cern.ch
//-----------------------------------------------------------------

#include <TH1.h>
#include <TTree.h>
#include <TClonesArray.h>

#include <AliITSVertexer3DTapan.h>
#include <AliITSRecPoint.h>
#include <AliITSgeomTGeo.h>
#include <AliESDVertex.h>

ClassImp(AliITSVertexer3DTapan)

void AliITSVertexer3DTapan::LoadClusters(TTree *cTree) {
  //--------------------------------------------------------------------
  //This function loads the SPD clusters
  //--------------------------------------------------------------------
   TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
   TBranch *branch=cTree->GetBranch("ITSRecPoints");
   branch->SetAddress(&clusters);

   Int_t nentr=cTree->GetEntries(),nc1=0,nc2=0;
   for (Int_t i=0; i<nentr; i++) {
       if (!cTree->GetEvent(i)) continue;
       //
       //   Below:
       //   "alpha" is the angle from the global X-axis to the
       //    local GEANT X'-axis  ( rot[0]=cos(alpha) and rot[1]=sin(alpha) )
       //    "phi" is the angle from the global X-axis to the
       //         local cluster X"-axis
       //

       Int_t    lay,lad,det; AliITSgeomTGeo::GetModuleId(i,lay,lad,det);

       if (lay>2) break;  //load the SPD clusters only

       Int_t ncl=clusters->GetEntriesFast();
       Float_t hPhi=0.;
       while (ncl--) {
          AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
	  Float_t pos[3];
	  c->GetGlobalXYZ(pos);
          if (lay==1) {
	    /*             fX1[nc1]= r*cp - c->GetY()*sp;
             fY1[nc1]= r*sp + c->GetY()*cp;
             fZ1[nc1]= c->GetZ(); */
	    fX1[nc1] = pos[0]; fY1[nc1] = pos[1]; fZ1[nc1] = pos[2];
	     CalculatePhi(fX1[nc1], fY1[nc1], hPhi);
	     fPhi1[nc1]= hPhi;
             nc1++;
          } else {
            /* fX2[nc2]= r*cp - c->GetY()*sp;
             fY2[nc2]= r*sp + c->GetY()*cp;
             fZ2[nc2]= c->GetZ(); */
	    fX2[nc2] = pos[0]; fY2[nc2] = pos[1]; fZ2[nc2] = pos[2];
	     CalculatePhi(fX2[nc2], fY2[nc2], hPhi);
	     fPhi2[nc2]= hPhi;
             nc2++;
          }
       }
   }
   ficlu1 = nc1; ficlu2 = nc2;
   AliInfo(Form("Number of clusters: %d (first layer) and %d (second layer)",ficlu1,ficlu2));
}

AliESDVertex *AliITSVertexer3DTapan::FindVertexForCurrentEvent(TTree *cTree) {
  //
  // This function reconstructs ....
  //
  //
  LoadClusters(cTree);

  Double_t pos[3], postemp[3];
  Double_t sigpos[3]={0.,0.,0.};
  Int_t ncontr, ncontrtemp;
  Float_t cuts[3];
  Int_t vtxstatus=0;

  //....
  pos[0] = 0.;   pos[1] = 0.;   pos[2] = 0.;
  cuts[0]=1.;  cuts[1]=1.;  cuts[2]=20.;
  CalculateVertex3d1(pos, cuts, ncontr);
  if(ncontr==0) {
    pos[0] = 9999.;   pos[1] = 9999.;   pos[2] = 9999.;
    vtxstatus = -1;
  }
  AliInfo(Form("1st step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus));

  if(vtxstatus == 0) {
    ncontrtemp = ncontr; postemp[0] = pos[0];   postemp[1] = pos[1];   postemp[2] = pos[2]; 
    cuts[0]=0.3;  cuts[1]=0.3;  cuts[2]=1.;
    CalculateVertex3d1(pos, cuts, ncontr);
    if(ncontr==0) {
      ncontr = ncontrtemp; pos[0] = postemp[0];   pos[1] = postemp[1];   pos[2] = postemp[2];
      vtxstatus = 2;
    }
    AliInfo(Form("2nd step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus));
  }

  if(vtxstatus == 0) {
    ncontrtemp = ncontr; postemp[0] = pos[0];   postemp[1] = pos[1];   postemp[2] = pos[2]; 
    cuts[0]=0.25;  cuts[1]=0.25;  cuts[2]=1.0;
    CalculateVertex3d2(pos, cuts, ncontr, sigpos);
    if(ncontr==0) {
      ncontr = ncontrtemp; pos[0] = postemp[0];   pos[1] = postemp[1];   pos[2] = postemp[2];
      vtxstatus = 3;
    }
    AliInfo(Form("3rd step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus));
  }

  if(vtxstatus == 0) {
    ncontrtemp = ncontr; postemp[0] = pos[0];   postemp[1] = pos[1];   postemp[2] = pos[2]; 
    cuts[0]=0.1;  cuts[1]=0.1;  cuts[2]=0.2;
    CalculateVertex3d2(pos, cuts, ncontr, sigpos);
    if(ncontr==0) {
      ncontr = ncontrtemp; pos[0] = postemp[0];   pos[1] = postemp[1];   pos[2] = postemp[2];
      vtxstatus = 4;
    }
    AliInfo(Form("4th step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus));
  }
  AliInfo(Form("Final step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus));

  Double_t covma[6]={0.,0.,0.,0.,0.,0.};
  covma[0]=sigpos[0];
  covma[2]=sigpos[1];
  covma[5]=sigpos[2];
  return new AliESDVertex(pos,covma,(Double_t)vtxstatus,ncontr,"AliITSVertexer3DTapan");

}


void AliITSVertexer3DTapan::CalculateVertex3d1(Double_t pos[3], Float_t cuts[3], Int_t &ncontr) {
  //
  // This function reconstructs first two steps of vertex
  //

  Double_t p1[4], p2[4], p3[4], p4[4];
  Double_t pa[3], pb[3];
  Double_t hphi1, hphi2, hphi3, hphi4;

  ncontr = 0;
  Float_t  phicut = 1.0;
  Double_t distance;  Float_t  distancecut = 1.0;
  Int_t    ibin=20;  Float_t  ilow=-1.; Float_t  ihigh=1.; 
  Int_t    ibinz=400; Float_t  ilowz=-20.; Float_t  ihighz=20.;

  TH1F *hx = new TH1F("hx","",  ibin, ilow, ihigh);
  TH1F *hy = new TH1F("hy","",  ibin, ilow, ihigh);
  TH1F *hz = new TH1F("hz","",  ibinz,ilowz,ihighz);
  
  for (Int_t ip1=0; ip1<ficlu1; ip1++) {
    // Two points on layer1: p1 and p3
    p1[0] = fX1[ip1];   p1[1] = fY1[ip1];    p1[2] = fZ1[ip1];   
    p3[0] = fX1[ip1+1]; p3[1] = fY1[ip1+1];  p3[2] = fZ1[ip1+1]; 
    hphi1 = fPhi1[ip1]; hphi3 = fPhi1[ip1+1];
    
    for (Int_t ip2=0; ip2<ficlu2; ip2++) {
      // Two points on layer 2: p2 and p4
      p2[0] = fX2[ip2];   p2[1] = fY2[ip2];     p2[2] = fZ2[ip2]; 
      p4[0] = fX2[ip2+1]; p4[1] = fY2[ip2+1];   p4[2] = fZ2[ip2+1]; 
      hphi2 = fPhi2[ip2]; hphi4 = fPhi2[ip2+1];

      // First line is formed by p1-p2 and second line by p3-p4
      // We find two points on each line which form the closest distance of the two lines
      // pa[0],pa[1],pa[2]: points on line 1 and pb[0],pb[1],pb[2]: points on line 2
      // Next: Consider x, y and z to be less than cuts[0], cuts[1] and cuts[2], respectively

      if(TMath::Abs(hphi1-hphi2)<phicut && TMath::Abs(hphi3-hphi4)<phicut){
	CalculateLine(p1, p2, p3, p4, pa, pb);

	if (pa[0]>pos[0]-cuts[0] && pa[0]<pos[0]+cuts[0] && pa[1]>pos[1]-cuts[1] && pa[1]<pos[1]+cuts[1] && pa[2]>pos[2]-cuts[2] && pa[2]<pos[2]+cuts[2]){
	  distance = (TMath::Sqrt(pow((pa[0]-pb[0]),2) + pow((pa[1]-pb[1]),2) + pow((pa[2]-pb[2]),2)));
	  if(distance<distancecut){
	    hx->Fill(pa[0]); 		hy->Fill(pa[1]);		hz->Fill(pa[2]);
	    hx->Fill(pb[0]); 		hy->Fill(pb[1]);		hz->Fill(pb[2]);
	    ncontr++;
	  }
	}
      }

      // Third line is formed by p1-p4 and fourth line by p3-p2
      // We find two points on each line which form the closest distance of the two lines
      // pa[0],pa[1],pa[2]: points on line 3 and pb[0],pb[1],pb[2]: points on line 4
      // Next: Consider x, y and z to be less than cuts[0], cuts[1] and cuts[2], respectively      
      if(TMath::Abs(hphi1-hphi4)<phicut && TMath::Abs(hphi3-hphi2)<phicut) {
	CalculateLine(p1, p4, p3, p2, pa, pb);
	if (pa[0]>pos[0]-cuts[0] && pa[0]<pos[0]+cuts[0] && pa[1]>pos[1]-cuts[1] && pa[1]<pos[1]+cuts[1]){
	  distance = (TMath::Sqrt(pow((pa[0]-pb[0]),2) + pow((pa[1]-pb[1]),2) + pow((pa[2]-pb[2]),2)));
	  if(distance<distancecut){
	    hx->Fill(pa[0]); 		hy->Fill(pa[1]);		hz->Fill(pa[2]);
	    hx->Fill(pb[0]); 		hy->Fill(pb[1]);		hz->Fill(pb[2]);
	    ncontr++;
	  }
	}
      }
    }
  }

  Int_t maxbinx = hx->GetMaximumBin(); 
  Int_t maxbiny = hy->GetMaximumBin();
  Int_t maxbinz = hz->GetMaximumBin();
  pos[0] = ilow  + ((ihigh-ilow)/ibin)*maxbinx;
  pos[1] = ilow  + ((ihigh-ilow)/ibin)*maxbiny;
  pos[2] = ilowz + ((ihighz-ilowz)/ibinz)*maxbinz;
  hx->Delete();
  hy->Delete();
  hz->Delete();
}

void AliITSVertexer3DTapan::CalculateVertex3d2(Double_t pos[3], Float_t cuts[3], Int_t &ncontr, Double_t sigpos[3]) {
  //
  // This function reconstructs second two steps of vertex
  //

  Double_t p1[4], p2[4], p3[4], p4[4];
  Double_t pa[3], pb[3];
  Double_t hphi1, hphi2, hphi3, hphi4;

  ncontr = 0;
  Float_t  phicut = 0.3;
  Double_t distance; Float_t  distancecut = 1.0;

  Double_t vertx  =0.;   Double_t verty  =0.;   Double_t vertz  =0.;	 
  Double_t vertx2 =0.;   Double_t verty2 =0.;   Double_t vertz2 =0.;	 

  for (Int_t ip1=0; ip1<ficlu1; ip1++) {
    // Two points on layer1: p1 and p3
    p1[0] = fX1[ip1];   p1[1] = fY1[ip1];    p1[2] = fZ1[ip1];   
    p3[0] = fX1[ip1+1]; p3[1] = fY1[ip1+1];  p3[2] = fZ1[ip1+1]; 
    hphi1 = fPhi1[ip1]; hphi3 = fPhi1[ip1+1];
    
    for (Int_t ip2=0; ip2<ficlu2; ip2++) {
      // Two points on layer 2: p2 and p4
      p2[0] = fX2[ip2];   p2[1] = fY2[ip2];     p2[2] = fZ2[ip2]; 
      p4[0] = fX2[ip2+1]; p4[1] = fY2[ip2+1];   p4[2] = fZ2[ip2+1]; 
      hphi2 = fPhi2[ip2]; hphi4 = fPhi2[ip2+1];

      // First line is formed by p1-p2 and second line by p3-p4
      // We find two points on each line which form the closest distance of the two lines
      // pa[0],pa[1],pa[2] are the points on line 1 and pb[0],pb[1],pb[2] are the points on line 2
      // Next: Consider x, y and z to be less than cuts[0], cuts[1] and cuts[2], respectively

      if(TMath::Abs(hphi1-hphi2)<phicut && TMath::Abs(hphi3-hphi4)<phicut){
	CalculateLine(p1, p2, p3, p4, pa, pb);

	// We consider the points where x, y and z points are less than xcut, ycut and zcut, respectively
	if (pa[0]>pos[0]-cuts[0] && pa[0]<pos[0]+cuts[0] && pa[1]>pos[1]-cuts[1] && pa[1]<pos[1]+cuts[1] && pa[2]>pos[2]-cuts[2] && pa[2]<pos[2]+cuts[2]){
	  distance = (TMath::Sqrt(pow((pa[0]-pb[0]),2) + pow((pa[1]-pb[1]),2) + pow((pa[2]-pb[2]),2)));
	  if(distance<distancecut){
	    ncontr++;
	    vertx  = vertx  + pa[0];       verty  = verty  + pa[1];       vertz  = vertz  + pa[2];
	    vertx2 = vertx2 + pa[0]*pa[0]; verty2 = verty2 + pa[1]*pa[1]; vertz2 = vertz2 + pa[2]*pa[2];
	    ncontr++;
	    vertx  = vertx  + pb[0];       verty  = verty  + pb[1];       vertz  = vertz  + pb[2];
	    vertx2 = vertx2 + pb[0]*pb[0]; verty2 = verty2 + pb[1]*pb[1]; vertz2 = vertz2 + pb[2]*pb[2];
	  }
	}
      }
      
      // Third line is formed by p1-p4 and fourth line by p3-p2
      // We find two points on each line which form the closest distance of the two lines
      // pa[0],pa[1],pa[2] are the points on line 3 and pb[0],pb[1],pb[2] are the points on line 4
      // Next: Consider x, y and z to be less than cuts[0], cuts[1] and cuts[2], respectively      
      if(TMath::Abs(hphi1-hphi4)<phicut && TMath::Abs(hphi3-hphi2)<phicut) {
	
	CalculateLine(p1, p4, p3, p2, pa, pb);
	if (pa[0]>pos[0]-cuts[0] && pa[0]<pos[0]+cuts[0] && pa[1]>pos[1]-cuts[1] && pa[1]<pos[1]+cuts[1] && pa[2]>pos[2]-cuts[2] && pa[2]<pos[2]+cuts[2]){
	  distance = (TMath::Sqrt(pow((pa[0]-pb[0]),2) + pow((pa[1]-pb[1]),2) + pow((pa[2]-pb[2]),2)));
	  if(distance<distancecut){
	    ncontr++;
	    vertx  = vertx  + pa[0];       verty  = verty  + pa[1];       vertz  = vertz  + pa[2];
	    vertx2 = vertx2 + pa[0]*pa[0]; verty2 = verty2 + pa[1]*pa[1]; vertz2 = vertz2 + pa[2]*pa[2];
	    ncontr++;
	    vertx  = vertx  + pb[0];       verty  = verty  + pb[1];       vertz  = vertz  + pb[2];
	    vertx2 = vertx2 + pb[0]*pb[0]; verty2 = verty2 + pb[1]*pb[1]; vertz2 = vertz2 + pb[2]*pb[2];
	  }
	}
      }
    }
  }

  if(ncontr>0){
    pos[0]  = vertx/ncontr;  pos[1] = verty/ncontr;   pos[2] = vertz/ncontr;
    vertx2 = vertx2/ncontr; verty2 = verty2/ncontr; vertz2 = vertz2/ncontr;
    sigpos[0] = TMath::Sqrt(vertx2 - pos[0]*pos[0]); 
    sigpos[1] = TMath::Sqrt(verty2 - pos[1]*pos[1]);
    sigpos[2] = TMath::Sqrt(vertz2 - pos[2]*pos[2]);
  }
  ncontr = ncontr/2;
}

void AliITSVertexer3DTapan::CalculatePhi(Float_t fx, Float_t fy, Float_t & phi)
{
  //calculates phi
  Float_t ybyx, phi1;
  const Float_t kradian = 180./3.141592654;
  
  if(fx==0.)
    {
      if(fy>0.) phi = 90.;
      if(fy<0.) phi = 270.;
    }
  if(fx != 0.)
    {
      ybyx = fy/fx;
      if(ybyx < 0) ybyx = - ybyx;
      phi1 = TMath::ATan(ybyx)*kradian;
      if(fx > 0 && fy > 0) phi = phi1;        // 1st Quadrant
      if(fx < 0 && fy > 0) phi = 180 - phi1;  // 2nd Quadrant
      if(fx < 0 && fy < 0) phi = 180 + phi1;  // 3rd Quadrant
      if(fx > 0 && fy < 0) phi = 360 - phi1;  // 4th Quadrant
      
    }
  phi = phi/kradian;
}

void AliITSVertexer3DTapan::CalculateLine(Double_t p1[4], Double_t p2[4], Double_t p3[4], Double_t p4[4], Double_t pa[3], Double_t pb[3]) const{
  //calculates line
  Double_t p13x, p13y, p13z;
  Double_t p21x, p21y, p21z;
  Double_t p43x, p43y, p43z;
  Double_t d1343, d4321, d1321, d4343, d2121;
  Double_t numer, denom;
  Double_t mua,   mub;
  mua = 0.; mub = 0.;
  
  p13x = p1[0] - p3[0];
  p13y = p1[1] - p3[1];
  p13z = p1[2] - p3[2];

  p21x = p2[0] - p1[0];
  p21y = p2[1] - p1[1];
  p21z = p2[2] - p1[2];

  p43x = p4[0] - p3[0];
  p43y = p4[1] - p3[1];
  p43z = p4[2] - p3[2];

  d1343 = p13x * p43x + p13y * p43y + p13z * p43z;
  d4321 = p43x * p21x + p43y * p21y + p43z * p21z;
  d1321 = p13x * p21x + p13y * p21y + p13z * p21z;
  d4343 = p43x * p43x + p43y * p43y + p43z * p43z;
  d2121 = p21x * p21x + p21y * p21y + p21z * p21z;

  denom = d2121 * d4343 - d4321 * d4321;
  numer = d1343 * d4321 - d1321 * d4343;

  if(denom>0) mua = numer / denom;
  if(d4343>0) mub = (d1343 + d4321 * (mua)) / d4343;

  pa[0] = p1[0] + mua * p21x;
  pa[1] = p1[1] + mua * p21y;
  pa[2] = p1[2] + mua * p21z;

  pb[0] = p3[0] + mub * p43x;
  pb[1] = p3[1] + mub * p43y;
  pb[2] = p3[2] + mub * p43z;
}
Commit	Line	Data
eb35e591	1	/**************************************************************************
	2	* Copyright(c) 2006-2008, 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	//-----------------------------------------------------------------
	17	// AliITSVertexer3DTapan class
	18	// This is a class for the 3d vertex finding
	19	// Origin: Tapan Nayak, VECC-CERN, Tapan.Nayak@cern.ch
	20	//-----------------------------------------------------------------
	21
	22	#include <TH1.h>
	23	#include <TTree.h>
	24	#include <TClonesArray.h>
	25
	26	#include <AliITSVertexer3DTapan.h>
	27	#include <AliITSRecPoint.h>
37add1d1	28	#include <AliITSgeomTGeo.h>
eb35e591	29	#include <AliESDVertex.h>
	30
	31	ClassImp(AliITSVertexer3DTapan)
	32
308c2f7c	33	void AliITSVertexer3DTapan::LoadClusters(TTree *cTree) {
eb35e591	34	//--------------------------------------------------------------------
	35	//This function loads the SPD clusters
	36	//--------------------------------------------------------------------
	37	TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
	38	TBranch *branch=cTree->GetBranch("ITSRecPoints");
	39	branch->SetAddress(&clusters);
	40
	41	Int_t nentr=cTree->GetEntries(),nc1=0,nc2=0;
	42	for (Int_t i=0; i<nentr; i++) {
	43	if (!cTree->GetEvent(i)) continue;
	44	//
	45	// Below:
	46	// "alpha" is the angle from the global X-axis to the
	47	// local GEANT X'-axis ( rot[0]=cos(alpha) and rot[1]=sin(alpha) )
	48	// "phi" is the angle from the global X-axis to the
	49	// local cluster X"-axis
	50	//
	51
37add1d1	52	Int_t lay,lad,det; AliITSgeomTGeo::GetModuleId(i,lay,lad,det);
eb35e591	53
	54	if (lay>2) break; //load the SPD clusters only
	55
eb35e591	56	Int_t ncl=clusters->GetEntriesFast();
225cff59	57	Float_t hPhi=0.;
eb35e591	58	while (ncl--) {
	59	AliITSRecPoint c=(AliITSRecPoint)clusters->UncheckedAt(ncl);
	60	Float_t pos[3];
	61	c->GetGlobalXYZ(pos);
	62	if (lay==1) {
	63	/* fX1[nc1]= rcp - c->GetY()sp;
	64	fY1[nc1]= rsp + c->GetY()cp;
	65	fZ1[nc1]= c->GetZ(); */
	66	fX1[nc1] = pos[0]; fY1[nc1] = pos[1]; fZ1[nc1] = pos[2];
	67	CalculatePhi(fX1[nc1], fY1[nc1], hPhi);
	68	fPhi1[nc1]= hPhi;
	69	nc1++;
	70	} else {
	71	/* fX2[nc2]= rcp - c->GetY()sp;
	72	fY2[nc2]= rsp + c->GetY()cp;
	73	fZ2[nc2]= c->GetZ(); */
	74	fX2[nc2] = pos[0]; fY2[nc2] = pos[1]; fZ2[nc2] = pos[2];
	75	CalculatePhi(fX2[nc2], fY2[nc2], hPhi);
	76	fPhi2[nc2]= hPhi;
	77	nc2++;
	78	}
	79	}
	80	}
	81	ficlu1 = nc1; ficlu2 = nc2;
	82	AliInfo(Form("Number of clusters: %d (first layer) and %d (second layer)",ficlu1,ficlu2));
eb35e591	83	}
eb35e591	84
308c2f7c	85	AliESDVertex AliITSVertexer3DTapan::FindVertexForCurrentEvent(TTree cTree) {
eb35e591	86	//
	87	// This function reconstructs ....
	88	//
	89	//
308c2f7c	90	LoadClusters(cTree);
eb35e591	91
225cff59	92	Double_t pos[3], postemp[3];
225cff59	93	Double_t sigpos[3]={0.,0.,0.};
eb35e591	94	Int_t ncontr, ncontrtemp;
	95	Float_t cuts[3];
	96	Int_t vtxstatus=0;
	97
	98	//....
	99	pos[0] = 0.; pos[1] = 0.; pos[2] = 0.;
	100	cuts[0]=1.; cuts[1]=1.; cuts[2]=20.;
	101	CalculateVertex3d1(pos, cuts, ncontr);
	102	if(ncontr==0) {
	103	pos[0] = 9999.; pos[1] = 9999.; pos[2] = 9999.;
	104	vtxstatus = -1;
	105	}
	106	AliInfo(Form("1st step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus));
	107
	108	if(vtxstatus == 0) {
	109	ncontrtemp = ncontr; postemp[0] = pos[0]; postemp[1] = pos[1]; postemp[2] = pos[2];
	110	cuts[0]=0.3; cuts[1]=0.3; cuts[2]=1.;
	111	CalculateVertex3d1(pos, cuts, ncontr);
	112	if(ncontr==0) {
	113	ncontr = ncontrtemp; pos[0] = postemp[0]; pos[1] = postemp[1]; pos[2] = postemp[2];
	114	vtxstatus = 2;
	115	}
	116	AliInfo(Form("2nd step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus));
	117	}
	118
	119	if(vtxstatus == 0) {
	120	ncontrtemp = ncontr; postemp[0] = pos[0]; postemp[1] = pos[1]; postemp[2] = pos[2];
	121	cuts[0]=0.25; cuts[1]=0.25; cuts[2]=1.0;
	122	CalculateVertex3d2(pos, cuts, ncontr, sigpos);
	123	if(ncontr==0) {
	124	ncontr = ncontrtemp; pos[0] = postemp[0]; pos[1] = postemp[1]; pos[2] = postemp[2];
	125	vtxstatus = 3;
	126	}
	127	AliInfo(Form("3rd step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus));
	128	}
	129
	130	if(vtxstatus == 0) {
	131	ncontrtemp = ncontr; postemp[0] = pos[0]; postemp[1] = pos[1]; postemp[2] = pos[2];
	132	cuts[0]=0.1; cuts[1]=0.1; cuts[2]=0.2;
	133	CalculateVertex3d2(pos, cuts, ncontr, sigpos);
	134	if(ncontr==0) {
	135	ncontr = ncontrtemp; pos[0] = postemp[0]; pos[1] = postemp[1]; pos[2] = postemp[2];
	136	vtxstatus = 4;
	137	}
	138	AliInfo(Form("4th step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus));
	139	}
	140	AliInfo(Form("Final step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus));
	141
225cff59	142	Double_t covma[6]={0.,0.,0.,0.,0.,0.};
	143	covma[0]=sigpos[0];
	144	covma[2]=sigpos[1];
	145	covma[5]=sigpos[2];
	146	return new AliESDVertex(pos,covma,(Double_t)vtxstatus,ncontr,"AliITSVertexer3DTapan");
308c2f7c	147
eb35e591	148	}
	149
	150
	151	void AliITSVertexer3DTapan::CalculateVertex3d1(Double_t pos[3], Float_t cuts[3], Int_t &ncontr) {
	152	//
	153	// This function reconstructs first two steps of vertex
	154	//
	155
	156	Double_t p1[4], p2[4], p3[4], p4[4];
	157	Double_t pa[3], pb[3];
	158	Double_t hphi1, hphi2, hphi3, hphi4;
	159
	160	ncontr = 0;
	161	Float_t phicut = 1.0;
	162	Double_t distance; Float_t distancecut = 1.0;
	163	Int_t ibin=20; Float_t ilow=-1.; Float_t ihigh=1.;
	164	Int_t ibinz=400; Float_t ilowz=-20.; Float_t ihighz=20.;
	165
	166	TH1F *hx = new TH1F("hx","", ibin, ilow, ihigh);
	167	TH1F *hy = new TH1F("hy","", ibin, ilow, ihigh);
	168	TH1F *hz = new TH1F("hz","", ibinz,ilowz,ihighz);
	169
	170	for (Int_t ip1=0; ip1<ficlu1; ip1++) {
	171	// Two points on layer1: p1 and p3
	172	p1[0] = fX1[ip1]; p1[1] = fY1[ip1]; p1[2] = fZ1[ip1];
	173	p3[0] = fX1[ip1+1]; p3[1] = fY1[ip1+1]; p3[2] = fZ1[ip1+1];
	174	hphi1 = fPhi1[ip1]; hphi3 = fPhi1[ip1+1];
	175
	176	for (Int_t ip2=0; ip2<ficlu2; ip2++) {
	177	// Two points on layer 2: p2 and p4
	178	p2[0] = fX2[ip2]; p2[1] = fY2[ip2]; p2[2] = fZ2[ip2];
	179	p4[0] = fX2[ip2+1]; p4[1] = fY2[ip2+1]; p4[2] = fZ2[ip2+1];
	180	hphi2 = fPhi2[ip2]; hphi4 = fPhi2[ip2+1];
	181
	182	// First line is formed by p1-p2 and second line by p3-p4
	183	// We find two points on each line which form the closest distance of the two lines
	184	// pa[0],pa[1],pa[2]: points on line 1 and pb[0],pb[1],pb[2]: points on line 2
	185	// Next: Consider x, y and z to be less than cuts[0], cuts[1] and cuts[2], respectively
	186
	187	if(TMath::Abs(hphi1-hphi2)<phicut && TMath::Abs(hphi3-hphi4)<phicut){
	188	CalculateLine(p1, p2, p3, p4, pa, pb);
	189
	190	if (pa[0]>pos[0]-cuts[0] && pa[0]<pos[0]+cuts[0] && pa[1]>pos[1]-cuts[1] && pa[1]<pos[1]+cuts[1] && pa[2]>pos[2]-cuts[2] && pa[2]<pos[2]+cuts[2]){
	191	distance = (TMath::Sqrt(pow((pa[0]-pb[0]),2) + pow((pa[1]-pb[1]),2) + pow((pa[2]-pb[2]),2)));
	192	if(distance<distancecut){
	193	hx->Fill(pa[0]); hy->Fill(pa[1]); hz->Fill(pa[2]);
	194	hx->Fill(pb[0]); hy->Fill(pb[1]); hz->Fill(pb[2]);
	195	ncontr++;
	196	}
	197	}
	198	}
	199
	200	// Third line is formed by p1-p4 and fourth line by p3-p2
	201	// We find two points on each line which form the closest distance of the two lines
	202	// pa[0],pa[1],pa[2]: points on line 3 and pb[0],pb[1],pb[2]: points on line 4
	203	// Next: Consider x, y and z to be less than cuts[0], cuts[1] and cuts[2], respectively
	204	if(TMath::Abs(hphi1-hphi4)<phicut && TMath::Abs(hphi3-hphi2)<phicut) {
	205	CalculateLine(p1, p4, p3, p2, pa, pb);
	206	if (pa[0]>pos[0]-cuts[0] && pa[0]<pos[0]+cuts[0] && pa[1]>pos[1]-cuts[1] && pa[1]<pos[1]+cuts[1]){
	207	distance = (TMath::Sqrt(pow((pa[0]-pb[0]),2) + pow((pa[1]-pb[1]),2) + pow((pa[2]-pb[2]),2)));
	208	if(distance<distancecut){
	209	hx->Fill(pa[0]); hy->Fill(pa[1]); hz->Fill(pa[2]);
	210	hx->Fill(pb[0]); hy->Fill(pb[1]); hz->Fill(pb[2]);
	211	ncontr++;
212	}
213	}
214	}
215	}
216	}
217
218	Int_t maxbinx = hx->GetMaximumBin();
219	Int_t maxbiny = hy->GetMaximumBin();
220	Int_t maxbinz = hz->GetMaximumBin();
221	pos[0] = ilow + ((ihigh-ilow)/ibin)*maxbinx;
222	pos[1] = ilow + ((ihigh-ilow)/ibin)*maxbiny;
223	pos[2] = ilowz + ((ihighz-ilowz)/ibinz)*maxbinz;
224	hx->Delete();
225	hy->Delete();
226	hz->Delete();
227	}
228
229	void AliITSVertexer3DTapan::CalculateVertex3d2(Double_t pos[3], Float_t cuts[3], Int_t &ncontr, Double_t sigpos[3]) {
230	//
231	// This function reconstructs second two steps of vertex
232	//
233
234	Double_t p1[4], p2[4], p3[4], p4[4];
235	Double_t pa[3], pb[3];
236	Double_t hphi1, hphi2, hphi3, hphi4;
237
238	ncontr = 0;
239	Float_t phicut = 0.3;
240	Double_t distance; Float_t distancecut = 1.0;
241
242	Double_t vertx =0.; Double_t verty =0.; Double_t vertz =0.;
243	Double_t vertx2 =0.; Double_t verty2 =0.; Double_t vertz2 =0.;
244
245	for (Int_t ip1=0; ip1<ficlu1; ip1++) {
246	// Two points on layer1: p1 and p3
247	p1[0] = fX1[ip1]; p1[1] = fY1[ip1]; p1[2] = fZ1[ip1];
248	p3[0] = fX1[ip1+1]; p3[1] = fY1[ip1+1]; p3[2] = fZ1[ip1+1];
249	hphi1 = fPhi1[ip1]; hphi3 = fPhi1[ip1+1];
250
251	for (Int_t ip2=0; ip2<ficlu2; ip2++) {
252	// Two points on layer 2: p2 and p4
253	p2[0] = fX2[ip2]; p2[1] = fY2[ip2]; p2[2] = fZ2[ip2];
254	p4[0] = fX2[ip2+1]; p4[1] = fY2[ip2+1]; p4[2] = fZ2[ip2+1];
255	hphi2 = fPhi2[ip2]; hphi4 = fPhi2[ip2+1];
256
257	// First line is formed by p1-p2 and second line by p3-p4
258	// We find two points on each line which form the closest distance of the two lines
259	// pa[0],pa[1],pa[2] are the points on line 1 and pb[0],pb[1],pb[2] are the points on line 2
260	// Next: Consider x, y and z to be less than cuts[0], cuts[1] and cuts[2], respectively
261
262	if(TMath::Abs(hphi1-hphi2)<phicut && TMath::Abs(hphi3-hphi4)<phicut){
263	CalculateLine(p1, p2, p3, p4, pa, pb);
264
265	// We consider the points where x, y and z points are less than xcut, ycut and zcut, respectively
266	if (pa[0]>pos[0]-cuts[0] && pa[0]<pos[0]+cuts[0] && pa[1]>pos[1]-cuts[1] && pa[1]<pos[1]+cuts[1] && pa[2]>pos[2]-cuts[2] && pa[2]<pos[2]+cuts[2]){
267	distance = (TMath::Sqrt(pow((pa[0]-pb[0]),2) + pow((pa[1]-pb[1]),2) + pow((pa[2]-pb[2]),2)));
268	if(distance<distancecut){
269	ncontr++;
270	vertx = vertx + pa[0]; verty = verty + pa[1]; vertz = vertz + pa[2];
271	vertx2 = vertx2 + pa[0]pa[0]; verty2 = verty2 + pa[1]pa[1]; vertz2 = vertz2 + pa[2]*pa[2];
272	ncontr++;
273	vertx = vertx + pb[0]; verty = verty + pb[1]; vertz = vertz + pb[2];
274	vertx2 = vertx2 + pb[0]pb[0]; verty2 = verty2 + pb[1]pb[1]; vertz2 = vertz2 + pb[2]*pb[2];
275	}
276	}
277	}
278
279	// Third line is formed by p1-p4 and fourth line by p3-p2
280	// We find two points on each line which form the closest distance of the two lines
281	// pa[0],pa[1],pa[2] are the points on line 3 and pb[0],pb[1],pb[2] are the points on line 4
282	// Next: Consider x, y and z to be less than cuts[0], cuts[1] and cuts[2], respectively
283	if(TMath::Abs(hphi1-hphi4)<phicut && TMath::Abs(hphi3-hphi2)<phicut) {
284
285	CalculateLine(p1, p4, p3, p2, pa, pb);
286	if (pa[0]>pos[0]-cuts[0] && pa[0]<pos[0]+cuts[0] && pa[1]>pos[1]-cuts[1] && pa[1]<pos[1]+cuts[1] && pa[2]>pos[2]-cuts[2] && pa[2]<pos[2]+cuts[2]){
287	distance = (TMath::Sqrt(pow((pa[0]-pb[0]),2) + pow((pa[1]-pb[1]),2) + pow((pa[2]-pb[2]),2)));
288	if(distance<distancecut){
289	ncontr++;
290	vertx = vertx + pa[0]; verty = verty + pa[1]; vertz = vertz + pa[2];
291	vertx2 = vertx2 + pa[0]pa[0]; verty2 = verty2 + pa[1]pa[1]; vertz2 = vertz2 + pa[2]*pa[2];
292	ncontr++;
293	vertx = vertx + pb[0]; verty = verty + pb[1]; vertz = vertz + pb[2];
294	vertx2 = vertx2 + pb[0]pb[0]; verty2 = verty2 + pb[1]pb[1]; vertz2 = vertz2 + pb[2]*pb[2];
295	}
296	}
297	}
298	}
299	}
300
301	if(ncontr>0){
302	pos[0] = vertx/ncontr; pos[1] = verty/ncontr; pos[2] = vertz/ncontr;
303	vertx2 = vertx2/ncontr; verty2 = verty2/ncontr; vertz2 = vertz2/ncontr;
304	sigpos[0] = TMath::Sqrt(vertx2 - pos[0]*pos[0]);
305	sigpos[1] = TMath::Sqrt(verty2 - pos[1]*pos[1]);
306	sigpos[2] = TMath::Sqrt(vertz2 - pos[2]*pos[2]);
307	}
308	ncontr = ncontr/2;
309	}
310
311	void AliITSVertexer3DTapan::CalculatePhi(Float_t fx, Float_t fy, Float_t & phi)
312	{
313	//calculates phi
314	Float_t ybyx, phi1;
315	const Float_t kradian = 180./3.141592654;
316
317	if(fx==0.)
318	{
319	if(fy>0.) phi = 90.;
320	if(fy<0.) phi = 270.;
321	}
322	if(fx != 0.)
323	{
324	ybyx = fy/fx;
325	if(ybyx < 0) ybyx = - ybyx;
326	phi1 = TMath::ATan(ybyx)*kradian;
327	if(fx > 0 && fy > 0) phi = phi1; // 1st Quadrant
328	if(fx < 0 && fy > 0) phi = 180 - phi1; // 2nd Quadrant
329	if(fx < 0 && fy < 0) phi = 180 + phi1; // 3rd Quadrant
330	if(fx > 0 && fy < 0) phi = 360 - phi1; // 4th Quadrant
331
332	}
333	phi = phi/kradian;
334	}
335
336	void AliITSVertexer3DTapan::CalculateLine(Double_t p1[4], Double_t p2[4], Double_t p3[4], Double_t p4[4], Double_t pa[3], Double_t pb[3]) const{
337	//calculates line
338	Double_t p13x, p13y, p13z;
339	Double_t p21x, p21y, p21z;
340	Double_t p43x, p43y, p43z;
341	Double_t d1343, d4321, d1321, d4343, d2121;
342	Double_t numer, denom;
343	Double_t mua, mub;
344	mua = 0.; mub = 0.;
345
346	p13x = p1[0] - p3[0];
347	p13y = p1[1] - p3[1];
348	p13z = p1[2] - p3[2];
349
350	p21x = p2[0] - p1[0];
351	p21y = p2[1] - p1[1];
352	p21z = p2[2] - p1[2];
353
354	p43x = p4[0] - p3[0];
355	p43y = p4[1] - p3[1];
356	p43z = p4[2] - p3[2];
357
358	d1343 = p13x * p43x + p13y * p43y + p13z * p43z;
359	d4321 = p43x * p21x + p43y * p21y + p43z * p21z;
360	d1321 = p13x * p21x + p13y * p21y + p13z * p21z;
361	d4343 = p43x * p43x + p43y * p43y + p43z * p43z;
362	d2121 = p21x * p21x + p21y * p21y + p21z * p21z;
363
364	denom = d2121 * d4343 - d4321 * d4321;
365	numer = d1343 * d4321 - d1321 * d4343;
366
367	if(denom>0) mua = numer / denom;
368	if(d4343>0) mub = (d1343 + d4321 * (mua)) / d4343;
369
370	pa[0] = p1[0] + mua * p21x;
371	pa[1] = p1[1] + mua * p21y;
372	pa[2] = p1[2] + mua * p21z;
373
374	pb[0] = p3[0] + mub * p43x;
375	pb[1] = p3[1] + mub * p43y;
376	pb[2] = p3[2] + mub * p43z;
377	}
378